library(readr)
library(taRifx)
library(tidyr)
library(plotly)
Loading required package: ggplot2
Registered S3 methods overwritten by 'ggplot2':
  method         from 
  [.quosures     rlang
  c.quosures     rlang
  print.quosures rlang
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

    last_plot

The following object is masked from ‘package:taRifx’:

    distinct

The following object is masked from ‘package:stats’:

    filter

The following object is masked from ‘package:graphics’:

    layout
FreeE_Propanediol.df <- read_table2("/home/aarcher/Downloads/FreeE_Propanediol.txt")
Parsed with column specification:
cols(
  `Axial-coord` = col_double(),
  `Radial-coord` = col_double(),
  freeE = col_double()
)
FreeE_Propanal.df <- read_table2("/home/aarcher/Downloads/FreeE_Propanal[2].txt")
Parsed with column specification:
cols(
  `Axial-coord` = col_double(),
  `Radial-coord` = col_double(),
  freeE = col_double()
)
FreeE_Propanediol.df <- FreeE_Propanediol.df[FreeE_Propanediol.df$`Axial-coord` > -9& FreeE_Propanediol.df$`Axial-coord` < 9,]
FreeE_Propanal.df <- FreeE_Propanal.df[FreeE_Propanal.df$`Axial-coord` > -9& FreeE_Propanal.df$`Axial-coord` < 9,]
FreeE_Propanediol_Center.df <- FreeE_Propanediol.df[FreeE_Propanediol.df$`Radial-coord` == 0.75 & FreeE_Propanediol.df$`Axial-coord` > -9& FreeE_Propanediol.df$`Axial-coord` < 9,]
FreeE_Propanal_Center.df <- FreeE_Propanal.df[FreeE_Propanal.df$`Radial-coord` == 0.75 & FreeE_Propanal.df$`Axial-coord` > -9& FreeE_Propanal.df$`Axial-coord` < 9,]
ggplot(data = FreeE_Propanediol.df, mapping = aes(x = `Radial-coord`, y = `Axial-coord`, fill = freeE)) + geom_tile() +  xlab(label = "Radius Coord")+  ylab(label = "Axial Coord") + ggtitle('Free Energy Distribution of 1,2-PDO')

plot(FreeE_Propanediol_Center.df$`Axial-coord`,FreeE_Propanediol_Center.df$freeE,ylab='Average Free Energy of 1,2-PDO', xlab='Axial Coord')

ggplot(data = FreeE_Propanal.df, mapping = aes(x = `Radial-coord`, y = `Axial-coord`, fill = freeE)) + geom_tile() +  xlab(label = "Radius Coord")+  ylab(label = "Axial Coord") + ggtitle('Free Energy Distribution of 1,2-PDO')

plot(FreeE_Propanal_Center.df$`Axial-coord`,FreeE_Propanal_Center.df$freeE,ylab='Average Free Energy of Proponal', xlab='Axial Coord')

corr <- function(tv1,tv2,eta1){
  K = matrix(1,nrow=dim(tv1)[1],ncol=dim(tv2)[1])
  for(dlcv in 1:dim(tv1)[2]){
  K <- K*exp(-eta1[dlcv]*outer(tv1[,dlcv],tv2[,dlcv],'-')^2)
  }
  K
}
MLE <-function(logeta, to, fo){
  m <- length(logeta)
  logeta1 <- logeta[-m]
  logeta2 <- logeta[m]
  
  # obtain the covariance via kronecker product
  
  corr0 <- corr(to, to, exp(logeta1))# c(., .)
  kappa0 <- corr0 + exp(logeta2)*diag(length(as.vector(fo)))
  kappa0_inv <- solve(kappa0)
  kappa0_inv_f0 <- solve(kappa0, as.vector(fo))
  kappa0_inv_1 <- solve(kappa0, rep(1, length(as.vector(fo))))# estimated values
  gammahat <- sum(kappa0_inv_f0)/sum(kappa0_inv_1)
  sigmahat <- mean(t(kappa0_inv_f0 - gammahat*kappa0_inv_1)*(as.vector(fo) - gammahat))
  return(list(gammahat=gammahat, sigmahat=sigmahat))
  }
# train each f as GP
negloglik <-function(logeta, t0, f0){
  n <- length(as.vector(f0))# obtain sigmahat for a given eta
  m <- length(logeta)
  logeta1 <- logeta[-m]
  logeta2 <- logeta[m]
  estimators <- MLE(logeta, t0, f0)
  sigmahat <- estimators$sigmahat# obtain the covariance matrices
  corr0 <- corr(t0, t0, exp(logeta1))# c(t0, t0)
  kappa0 <- corr0 + exp(logeta2)*diag(n)
  # obtain the log of determinants efficiently
  logdetkappa0 <- determinant(kappa0,logarithm=TRUE)$modulus
  return(1/2*logdetkappa0 + n/2*log(sigmahat))
  }
fitGP <-function(logeta=c(1, 1, 1), lowerb=c(-1, -1,-1), upperb=c(4, 4,4), t_tr, f_tr){
  log_eta_hat <- optim(par=logeta,fn=negloglik,lower=lowerb,upper=upperb,method="L-BFGS-B", t0=t_tr, f0=f_tr)$par
  estimators <- MLE(log_eta_hat, t_tr, f_tr)
  m <- length(log_eta_hat)
  return(list(etahat=exp(log_eta_hat[-m]),g=exp(log_eta_hat[m]),sigmahat=estimators$sigmahat,gammahat=estimators$gammahat))
}
predictGP <-function(fitted_info, t_test, t_tr, f_tr){# obtain the fitted info
  etahat1 <- fitted_info$etahat
  g <- fitted_info$g
  sigmahat <- fitted_info$sigmahat
  gammahat <- fitted_info$gammahat
  
  # obtain the final covariance structure
  
  corr0 <- corr(t_tr, t_tr, etahat1)  # c(t_tr, t_tr)
  kappa0 <- corr0 + g*diag(length(as.vector(f_tr)))
  kappa0_inv <- (1/sigmahat)*solve(kappa0)
  kappa0_inv_f0 <- kappa0_inv%*%(as.vector(f_tr) - gammahat)
  corr_Tv <- corr(t_test, t_tr, etahat1)# c(t_test, t_tr)
  kappa_Tv <- sigmahat*corr_Tv# c(t_test, t_tr)
  
  corr_toTv <- corr(t_tr, t_test, etahat1)# c(t_tr, t_test)
  kappa_toTv <- sigmahat *corr_toTv#  c(., .)
  
  gamma_vec <- (gammahat + kappa_Tv%*%kappa0_inv_f0)
  dkappa <- sigmahat * rep(1, dim(t_test)[1]*dim(f_tr)[2])
  kappa_vec <- (pmax(dkappa -apply(t((kappa_Tv)%*%(kappa0_inv))*(kappa_toTv), 2, sum), 10**(-13)))
  return(list(pred_mean=gamma_vec, pred_var=kappa_vec))}
ftrain <- as.matrix(FreeE_Propanal.df[,3])
inputtrain <-as.matrix(FreeE_Propanal.df[,1:2])
fitted_info <- fitGP(logeta=c(1, 1,1), lowerb=c(0, 0,0), upperb=c(5, 5,5), inputtrain, ftrain)
fitted_info
$etahat
[1] 1 1

$g
[1] 1

$sigmahat
[1] 2.181122

$gammahat
[1] 5.102003
predmean =rep(0,dim(inputtrain)[1])
predLB =  rep(0,dim(inputtrain)[1])
predUB =rep(0,dim(inputtrain)[1])
for(k in 1:dim(inputtrain)[1]){
  predinfo <- predictGP(fitted_info,t(as.matrix(inputtrain[k,])), inputtrain, ftrain)
  predmean[k] <- predinfo$pred_mean
  predLB[k] = predmean[k]- qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUB[k] = predmean[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(paste('PCGPMethod: Test MSE:',round(mean((predmean-ftrain)^2),3),', Coverage:',100*round(mean((predLB<ftrain)*(predUB>ftrain)),3),'%',sep =''))
[1] "PCGPMethod: Test MSE:0.764, Coverage:95.6%"
FreeE_Propanal_Pred.df <- FreeE_Propanal.df
FreeE_Propanal_LB.df <- FreeE_Propanal.df
FreeE_Propanal_UB.df <- FreeE_Propanal.df
FreeE_Propanal_Pred.df$freeE <- predmean
FreeE_Propanal_LB.df$freeE <- predLB
FreeE_Propanal_UB.df$freeE <- predUB
z <- as.matrix(spread(FreeE_Propanal.df, `Radial-coord`, freeE)[,-1])
zpred <- as.matrix(spread(FreeE_Propanal_Pred.df, `Radial-coord`, freeE)[,-1])
zLB <- as.matrix(spread(FreeE_Propanal_LB.df, `Radial-coord`, freeE)[,-1])
zUB <- as.matrix(spread(FreeE_Propanal_UB.df, `Radial-coord`, freeE)[,-1])
color <- rep(1, length(z))
dim(color) <- dim(z)
x <-unique(FreeE_Propanal.df$`Radial-coord`)
y <-unique(FreeE_Propanal.df$`Axial-coord`)
p<-  plot_ly() %>% add_trace(x=FreeE_Propanal.df$`Radial-coord`, y=FreeE_Propanal.df$`Axial-coord`, z=FreeE_Propanal.df$freeE, mode= "markers",type = "scatter3d",marker = list(size = 2, color = "blue", symbol = 104), opacity=0.7,showlegend = FALSE )
p <- p %>% add_trace( x=x, y=y, z=zpred,type = "surface",colorscale = list(c(0, 1), c("tan", "blue")),opacity=0.7)  
df1 <- split(FreeE_Propanal_LB.df, FreeE_Propanal_LB.df$`Radial-coord`)
df2 <- split(FreeE_Propanal_LB.df, FreeE_Propanal_LB.df$`Axial-coord`)
for(i in seq_along(df1)){
  df_sp <- df1[[i]]
  p <- add_trace(p, line = list(
    color = "#000000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
for(i in seq_along(df2)){
  df_sp <- df2[[i]]
  p <- add_trace(p, line = list(
    color = "#000000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
df1 <- split(FreeE_Propanal_UB.df, FreeE_Propanal_UB.df$`Radial-coord`)
df2 <- split(FreeE_Propanal_UB.df, FreeE_Propanal_UB.df$`Axial-coord`)
for(i in seq_along(df1)){
  df_sp <- df1[[i]]
  p <- add_trace(p, line = list(
    color = "#FF0000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
for(i in seq_along(df2)){
  df_sp <- df2[[i]]
  p <- add_trace(p, line = list(
    color = "#FF0000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
p%>% layout(scene = list( xaxis = list(title = "Radius, r"),yaxis = list(title = "Axial, z")))
axiallen <- 1000
axialcoord <- unique(FreeE_Propanal.df$`Axial-coord`)
centercoords <-as.matrix( data.frame(seq(min(axialcoord),max(axialcoord),length.out=axiallen),rep(0,axiallen)))
predmean =rep(0,dim(centercoords)[1])
predLB =  rep(0,dim(centercoords)[1])
predUB =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGP(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmean[k] <- predinfo$pred_mean
  predLB[k] = predmean[k]- qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUB[k] = predmean[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
plot(centercoords[,1],predmean, ylim = c(min(predLB),max(predUB)))
lines(centercoords[,1],predLB)
lines(centercoords[,1],predUB)

deltax = (max(axialcoord)-min(axialcoord))/(axiallen-1)
predmeandev = (predmean[3:axiallen]-predmean[1:(axiallen-2)])/deltax
preddevLB = (predLB[3:axiallen]-predLB[1:(axiallen-2)])/deltax
preddevUB = (predUB[3:axiallen]-predUB[1:(axiallen-2)])/deltax
predmeandevdev = (predmean[3:axiallen]-2*predmean[2:(axiallen-1)] + predmean[1:(axiallen-2)])/deltax^2
preddevdevLB = (predLB[3:axiallen]-2*predmean[2:(axiallen-1)] +predLB[1:(axiallen-2)])/deltax^2
preddevdevUB = (predUB[3:axiallen]-2*predmean[2:(axiallen-1)] +predUB[1:(axiallen-2)])/deltax^2
indmean = which(predmean==max(predmean))
indLB = which(predLB==max(predLB))
indUB = which(predUB==max(predUB))
print(predmean[indmean])
[1] 5.922285
print(predLB[indLB])
[1] 3.288251
print(predUB[indUB])
[1] 8.559228
plot(centercoords[2:(axiallen-1),1],predmeandev , ylim = c(min(predmeandev,preddevLB,preddevUB),max(predmeandev,preddevLB,preddevUB)))
lines(centercoords[2:(axiallen-1),1],preddevLB)
lines(centercoords[2:(axiallen-1),1],preddevUB)

print(predmeandev[indmean])
[1] -0.01772556
print(preddevLB[indLB])
[1] -0.0243247
print(preddevUB[indUB])
[1] -0.01114983
plot(centercoords[2:(axiallen-1),1],predmeandevdev , ylim = c(min(predmeandevdev,preddevdevLB,preddevdevUB),max(predmeandevdev,preddevdevLB,preddevdevUB)))
lines(centercoords[2:(axiallen-1),1],preddevdevLB)
lines(centercoords[2:(axiallen-1),1],preddevdevUB)

print(predmeandevdev[indmean])
[1] -0.4683776
print(preddevdevLB[indLB])
[1] -18189.59
print(preddevdevUB[indUB])
[1] 18223.72
corrdevobv <- function(tv1,tv2,eta1){
  K <- -2*eta1[1]*exp(-eta1[1]*outer(tv1[,1],tv2[,1],'-')^2)*(outer(tv1[,1],tv2[,1],'-'))
  K <- K*exp(-eta1[2]*abs(outer(tv1[,2],tv2[,2],'-'))^2)
  K
}
predictGPdev <-function(fitted_info, t_test, t_tr, f_tr){# obtain the fitted info
  etahat1 <- fitted_info$etahat
  g <- fitted_info$g
  sigmahat <- fitted_info$sigmahat
  gammahat <- fitted_info$gammahat
  
  # obtain the final covariance structure
  
  corr0 <- corr(t_tr, t_tr, etahat1)  # c(t_tr, t_tr)
  kappa0 <- corr0 + g*diag(length(as.vector(f_tr)))
  kappa0_inv <- (1/sigmahat)*solve(kappa0)
  kappa0_inv_f0 <- kappa0_inv%*%(as.vector(f_tr) - gammahat)
  corr_Tv <- corrdevobv(t_test, t_tr, etahat1)# c(t_test, t_tr)
  kappa_Tv <- sigmahat*corr_Tv# c(t_test, t_tr)
  
  corr_toTv <- corrdevobv(t_tr, t_test, etahat1)# c(t_tr, t_test)
  kappa_toTv <- sigmahat *corr_toTv#  c(., .)
  
  gamma_vec <-  kappa_Tv%*%kappa0_inv_f0
  dkappa <- sigmahat*(1 + g)*rep(1, dim(t_test)[1]*dim(f_tr)[2])
  kappa_vec <- (pmax(dkappa -apply(t((kappa_Tv)%*%(kappa0_inv))*(kappa_toTv), 2, sum), 10**(-13)))
  return(list(pred_mean=gamma_vec, pred_var=kappa_vec))}
predmeandev =rep(0,dim(centercoords)[1])
predLBdev =  rep(0,dim(centercoords)[1])
predUBdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandev[k] <- predinfo$pred_mean
  predLBdev[k] = predmeandev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdev[k] = predmeandev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
plot(centercoords[,1],predmeandev)# ylim = c(min(predLBdev),max(predUBdev)))
lines(centercoords[,1],predLBdev)
lines(centercoords[,1],predUBdev)

inds = which(diff(sign(predmeandev)) != 0)
print(predmeandev[inds])
[1] -0.0013450402  0.0007428134 -0.0036200835  0.0072476454 -0.0017144610
print(predLBdev[inds])
[1] -4.333562 -4.313164 -4.325467 -4.294252 -4.301100
print(predUBdev[inds])
[1] 4.330872 4.314650 4.318227 4.308747 4.297671
corrdevdevobv <- function(tv1,tv2,eta1){
  K <- 2*eta1[1]*exp(-eta1[1]*outer(tv1[,1],tv2[,1],'-')^2)*(2*eta1[1]*(outer(tv1[,1],tv2[,1],'-'))^2-1)
  K <- K*exp(-eta1[2]*abs(outer(tv1[,2],tv2[,2],'-'))^2)
  K
}
corrdevdevdevobv <- function(tv1,tv2,eta1){
  K <- -4*(eta1[1]^2)*(outer(tv1[,1],tv2[,1],'-'))*exp(-eta1[1]*outer(tv1[,1],tv2[,1],'-')^2)*(2*eta1[1]*(outer(tv1[,1],tv2[,1],'-'))^2-3)
  K <- K*exp(-eta1[2]*abs(outer(tv1[,2],tv2[,2],'-'))^2)
  K
}
corrdevdevdevdevobv <- function(tv1,tv2,eta1){
  K <- 4*(eta1[1]^2)*exp(-eta1[1]*outer(tv1[,1],tv2[,1],'-')^2)*(4*eta1[1]*(outer(tv1[,1],tv2[,1],'-')^2)*(eta1[1]*(outer(tv1[,1],tv2[,1],'-'))^2-3)+3)
  K <- K*exp(-eta1[2]*abs(outer(tv1[,2],tv2[,2],'-'))^2)
  K
}
predictGPdevdev <-function(fitted_info, t_test, t_tr, f_tr){# obtain the fitted info
  etahat1 <- fitted_info$etahat
  g <- fitted_info$g
  sigmahat <- fitted_info$sigmahat
  gammahat <- fitted_info$gammahat
  
  # obtain the final covariance structure
  
  corr0 <- corr(t_tr, t_tr, etahat1)  # c(t_tr, t_tr)
  kappa0 <- corr0 + g*diag(length(as.vector(f_tr)))
  kappa0_inv <- (1/sigmahat)*solve(kappa0)
  kappa0_inv_f0 <- kappa0_inv%*%(as.vector(f_tr) - gammahat)
  corr_Tv <- corrdevdevobv(t_test, t_tr, etahat1)# c(t_test, t_tr)
  kappa_Tv <- sigmahat*corr_Tv# c(t_test, t_tr)
  
  corr_toTv <- corrdevdevobv(t_tr, t_test, etahat1)# c(t_tr, t_test)
  kappa_toTv <- sigmahat *corr_toTv#  c(., .)
  
  gamma_vec <-  kappa_Tv%*%kappa0_inv_f0
  dkappa <- sigmahat*(1 + g)*rep(1, dim(t_test)[1]*dim(f_tr)[2])
  kappa_vec <- (pmax(dkappa -apply(t((kappa_Tv)%*%(kappa0_inv))*(kappa_toTv), 2, sum), 10**(-13)))
  return(list(pred_mean=gamma_vec, pred_var=kappa_vec))}
predictGPdevdevdev <-function(fitted_info, t_test, t_tr, f_tr){# obtain the fitted info
  etahat1 <- fitted_info$etahat
  g <- fitted_info$g
  sigmahat <- fitted_info$sigmahat
  gammahat <- fitted_info$gammahat
  
  # obtain the final covariance structure
  
  corr0 <- corr(t_tr, t_tr, etahat1)  # c(t_tr, t_tr)
  kappa0 <- corr0 + g*diag(length(as.vector(f_tr)))
  kappa0_inv <- (1/sigmahat)*solve(kappa0)
  kappa0_inv_f0 <- kappa0_inv%*%(as.vector(f_tr) - gammahat)
  corr_Tv <- corrdevdevdevobv(t_test, t_tr, etahat1)# c(t_test, t_tr)
  kappa_Tv <- sigmahat*corr_Tv# c(t_test, t_tr)
  
  corr_toTv <- corrdevdevdevobv(t_tr, t_test, etahat1)# c(t_tr, t_test)
  kappa_toTv <- sigmahat *corr_toTv#  c(., .)
  
  gamma_vec <-  kappa_Tv%*%kappa0_inv_f0
  dkappa <- sigmahat*(1 + g)*rep(1, dim(t_test)[1]*dim(f_tr)[2])
  kappa_vec <- (pmax(dkappa -apply(t((kappa_Tv)%*%(kappa0_inv))*(kappa_toTv), 2, sum), 10**(-13)))
  return(list(pred_mean=gamma_vec, pred_var=kappa_vec))}
predictGPdevdevdevdev <-function(fitted_info, t_test, t_tr, f_tr){# obtain the fitted info
  etahat1 <- fitted_info$etahat
  g <- fitted_info$g
  sigmahat <- fitted_info$sigmahat
  gammahat <- fitted_info$gammahat
  
  # obtain the final covariance structure
  
  corr0 <- corr(t_tr, t_tr, etahat1)  # c(t_tr, t_tr)
  kappa0 <- corr0 + g*diag(length(as.vector(f_tr)))
  kappa0_inv <- (1/sigmahat)*solve(kappa0)
  kappa0_inv_f0 <- kappa0_inv%*%(as.vector(f_tr) - gammahat)
  corr_Tv <- corrdevdevdevdevobv(t_test, t_tr, etahat1)# c(t_test, t_tr)
  kappa_Tv <- sigmahat*corr_Tv# c(t_test, t_tr)
  
  corr_toTv <- corrdevdevdevdevobv(t_tr, t_test, etahat1)# c(t_tr, t_test)
  kappa_toTv <- sigmahat *corr_toTv#  c(., .)
  
  gamma_vec <-  kappa_Tv%*%kappa0_inv_f0
  dkappa <- sigmahat*(1 + g)*rep(1, dim(t_test)[1]*dim(f_tr)[2])
  kappa_vec <- (pmax(dkappa -apply(t((kappa_Tv)%*%(kappa0_inv))*(kappa_toTv), 2, sum), 10**(-13)))
  return(list(pred_mean=gamma_vec, pred_var=kappa_vec))}
predmeandevdev =rep(0,dim(centercoords)[1])
predLBdevdev =  rep(0,dim(centercoords)[1])
predUBdevdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdevdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandevdev[k] <- predinfo$pred_mean
  predLBdevdev[k] = predmeandevdev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdevdev[k] = predmeandevdev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(predmeandevdev[inds])
[1]  0.4255432 -0.1365788  0.3426273 -0.4790840  0.1868095
print(predLBdevdev[inds])
[1] -3.365105 -3.459251 -3.072764 -3.648487 -2.956790
print(predUBdevdev[inds])
[1] 4.216191 3.186093 3.758018 2.690319 3.330409
plot(centercoords[,1],predmeandevdev, ylim = c(min(predLBdev),max(predUBdev)))
lines(centercoords[,1],predLBdevdev)
lines(centercoords[,1],predUBdevdev)

predmeandevdevdev =rep(0,dim(centercoords)[1])
predLBdevdevdev =  rep(0,dim(centercoords)[1])
predUBdevdevdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdevdevdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandevdevdev[k] <- predinfo$pred_mean
  predLBdevdevdev[k] = predmeandevdevdev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdevdevdev[k] = predmeandevdevdev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(predmeandevdevdev[inds])
[1] -0.95180132  0.09850211  0.53404848  0.29814148 -0.26788548
print(predLBdevdevdev[inds])
[1] -9.368393 -7.376385 -7.329279 -6.538623 -6.984487
print(predUBdevdevdev[inds])
[1] 7.464790 7.573389 8.397376 7.134906 6.448717
plot(centercoords[,1],predmeandevdevdev, ylim = c(min(predLBdevdevdev),max(predUBdevdevdev)))
lines(centercoords[,1],predLBdevdevdev)
lines(centercoords[,1],predUBdevdevdev)

predmeandevdevdevdev =rep(0,dim(centercoords)[1])
predLBdevdevdevdev =  rep(0,dim(centercoords)[1])
predUBdevdevdevdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdevdevdevdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandevdevdevdev[k] <- predinfo$pred_mean
  predLBdevdevdevdev[k] = predmeandevdevdevdev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdevdevdevdev[k] = predmeandevdevdevdev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(predmeandevdevdevdev[inds])
[1] -0.8613916  0.2039218 -0.8864295  0.9720034 -0.2326909
print(predLBdevdevdevdev[inds])
[1] -2.5334432  0.2039211 -0.8864301  0.9720028 -0.2326915
print(predUBdevdevdevdev[inds])
[1]  0.8106600  0.2039224 -0.8864289  0.9720040 -0.2326903
plot(centercoords[,1],predmeandevdevdevdev, ylim = c(min(predLBdevdevdevdev),max(predUBdevdevdevdev)))
lines(centercoords[,1],predLBdevdevdevdev)
lines(centercoords[,1],predUBdevdevdevdev)

kBT = (293*0.001985)
print((exp(predmean[inds[4]]-min(predmean))/kBT)*sqrt(-2*pi*kBT/predmeandevdev[inds[4]]) + kBT*((1/8)*(predmeandevdevdevdev[inds[4]]/(predmeandevdev[inds[4]])^2) - (5/24)*((predmeandevdevdev[inds[4]])^2/(predmeandevdev[inds[4]])^3)))
[1] 24.22749
print((exp(predLB[inds[4]]-min(predLB))/kBT)*sqrt(-2*pi*kBT/predLBdevdev[inds[4]])+ kBT*((1/8)*(predLBdevdevdevdev[inds[4]]/(predLBdevdev[inds[4]])^2) - (5/24)*((predLBdevdevdev[inds[4]])^2/(predLBdevdev[inds[4]])^3)))
[1] 8.784958
ftrain <- as.matrix(FreeE_Propanediol.df[,3])
inputtrain <-as.matrix(FreeE_Propanediol.df[,1:2])
fitted_info <- fitGP(logeta=c(1, 1,1), lowerb=c(0, 0,0), upperb=c(5, 5,5), inputtrain, ftrain)
fitted_info
$etahat
[1] 1 1

$g
[1] 1

$sigmahat
[1] 3.769055

$gammahat
[1] 6.751834
predmean =rep(0,dim(inputtrain)[1])
predLB =  rep(0,dim(inputtrain)[1])
predUB =rep(0,dim(inputtrain)[1])
for(k in 1:dim(inputtrain)[1]){
  predinfo <- predictGP(fitted_info,t(as.matrix(inputtrain[k,])), inputtrain, ftrain)
  predmean[k] <- predinfo$pred_mean
  predLB[k] = predmean[k]- qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUB[k] = predmean[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(paste('PCGPMethod: Test MSE:',round(mean((predmean-ftrain)^2),3),', Coverage:',100*round(mean((predLB<ftrain)*(predUB>ftrain)),3),'%',sep =''))
[1] "PCGPMethod: Test MSE:1.294, Coverage:100%"
FreeE_Propanediol_Pred.df <- FreeE_Propanediol.df
FreeE_Propanediol_LB.df <- FreeE_Propanediol.df
FreeE_Propanediol_UB.df <- FreeE_Propanediol.df
FreeE_Propanediol_Pred.df$freeE <- predmean
FreeE_Propanediol_LB.df$freeE <- predLB
FreeE_Propanediol_UB.df$freeE <- predUB
z <- as.matrix(spread(FreeE_Propanediol.df, `Radial-coord`, freeE)[,-1])
zpred <- as.matrix(spread(FreeE_Propanediol_Pred.df, `Radial-coord`, freeE)[,-1])
zLB <- as.matrix(spread(FreeE_Propanediol_LB.df, `Radial-coord`, freeE)[,-1])
zUB <- as.matrix(spread(FreeE_Propanediol_UB.df, `Radial-coord`, freeE)[,-1])
color <- rep(1, length(z))
dim(color) <- dim(z)
x <-unique(FreeE_Propanediol.df$`Radial-coord`)
y <-unique(FreeE_Propanediol.df$`Axial-coord`)
p <- plot_ly() %>% add_trace( x=x, y=y, z=z,type = "surface",colorscale = list(c(0, 1), c("tan", "blue")))
df1 <- split(FreeE_Propanediol_LB.df, FreeE_Propanediol_LB.df$`Radial-coord`)
df2 <- split(FreeE_Propanediol_LB.df, FreeE_Propanediol_LB.df$`Axial-coord`)
for(i in seq_along(df1)){
  df_sp <- df1[[i]]
  p <- add_trace(p, line = list(
    color = "#000000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
for(i in seq_along(df2)){
  df_sp <- df2[[i]]
  p <- add_trace(p, line = list(
    color = "#000000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
df1 <- split(FreeE_Propanediol_UB.df, FreeE_Propanediol_UB.df$`Radial-coord`)
df2 <- split(FreeE_Propanediol_UB.df, FreeE_Propanediol_UB.df$`Axial-coord`)
for(i in seq_along(df1)){
  df_sp <- df1[[i]]
  p <- add_trace(p, line = list(
    color = "#FF0000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
for(i in seq_along(df2)){
  df_sp <- df2[[i]]
  p <- add_trace(p, line = list(
    color = "#FF0000", 
    width = 2
  ), 
  mode = "lines", 
  type = "scatter3d", 
  x = df_sp$`Radial-coord`,
  y = df_sp$`Axial-coord`,
  z = df_sp$freeE,
  showlegend = FALSE)
}
p%>% layout(scene = list( xaxis = list(title = "Radius, r"),yaxis = list(title = "Axial, z")))
axiallen <- 1000
axialcoord <- unique(FreeE_Propanal.df$`Axial-coord`)
centercoords <-as.matrix( data.frame(seq(min(axialcoord),max(axialcoord),length.out=axiallen),rep(0,axiallen)))
predmean =rep(0,dim(centercoords)[1])
predLB =  rep(0,dim(centercoords)[1])
predUB =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGP(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmean[k] <- predinfo$pred_mean
  predLB[k] = predmean[k]- qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUB[k] = predmean[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
ind = which(predmean==max(predmean))
print(centercoords[ind,1])
seq.min.axialcoord...max.axialcoord...length.out...axiallen. 
                                                  -0.6891892 
print(predmean[ind])
[1] 6.308108
print(predLB[ind])
[1] 2.846148
print(predUB[ind])
[1] 9.770069
plot(centercoords[,1],predmean, ylim = c(min(predLB),max(predUB)))
lines(centercoords[,1],predLB)
lines(centercoords[,1],predUB)

predmeandev =rep(0,dim(centercoords)[1])
predLBdev =  rep(0,dim(centercoords)[1])
predUBdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandev[k] <- predinfo$pred_mean
  predLBdev[k] = predmeandev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdev[k] = predmeandev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
inds = which(diff(sign(predmeandev)) != 0)
print(centercoords[inds,1])
[1] -8.0405405 -7.5810811 -5.9134134 -0.6891892  2.1016016  4.4329329  6.9344344
print(predmeandev[inds])
[1] -0.0007317956  0.0033212283 -0.0024386434  0.0037187705 -0.0117369049  0.0007106318 -0.0009734004
print(predLBdev[inds])
[1] -5.698458 -5.618673 -5.702873 -5.641543 -5.709823 -5.620143 -5.704034
print(predUBdev[inds])
[1] 5.696995 5.625315 5.697995 5.648981 5.686349 5.621564 5.702087
plot(centercoords[,1],predmeandev, ylim = c(min(predLBdev),max(predUBdev)))
lines(centercoords[,1],predLBdev)
lines(centercoords[,1],predUBdev)

predmeandevdev =rep(0,dim(centercoords)[1])
predLBdevdev =  rep(0,dim(centercoords)[1])
predUBdevdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdevdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandevdev[k] <- predinfo$pred_mean
  predLBdevdev[k] = predmeandevdev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdevdev[k] = predmeandevdev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(predmeandevdev[inds])
[1]  0.3513365 -0.2631869  0.1510140 -0.5502032  0.7288531 -0.3181518  0.1565459
print(predLBdevdev[inds])
[1] -4.633672 -4.018359 -4.557604 -4.598919 -3.953395 -4.040190 -4.585034
print(predUBdevdev[inds])
[1] 5.336345 3.491985 4.859632 3.498513 5.411101 3.403887 4.898126
plot(centercoords[,1],predmeandevdev, ylim = c(min(predLBdev),max(predUBdev)))
lines(centercoords[,1],predLBdevdev)
lines(centercoords[,1],predUBdevdev)

predmeandevdevdev =rep(0,dim(centercoords)[1])
predLBdevdevdev =  rep(0,dim(centercoords)[1])
predUBdevdevdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdevdevdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandevdevdev[k] <- predinfo$pred_mean
  predLBdevdevdev[k] = predmeandevdevdev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdevdevdev[k] = predmeandevdevdev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(predmeandevdevdev[inds])
[1] -1.6050450 -0.7300791 -0.3220568 -0.4499744  0.4405389  0.3292321 -0.3360387
print(predLBdevdevdev[inds])
[1] -12.716097  -7.507184 -11.529364  -8.921627 -10.663563  -6.607474 -11.658531
print(predUBdevdevdev[inds])
[1]  9.506007  6.047025 10.885250  8.021678 11.544641  7.265938 10.986453
plot(centercoords[,1],predmeandevdevdev, ylim = c(min(predLBdevdevdev),max(predUBdevdevdev)))
lines(centercoords[,1],predLBdevdevdev)
lines(centercoords[,1],predUBdevdevdev)

predmeandevdevdevdev =rep(0,dim(centercoords)[1])
predLBdevdevdevdev =  rep(0,dim(centercoords)[1])
predUBdevdevdevdev =rep(0,dim(centercoords)[1])
for(k in 1:dim(centercoords)[1]){
  predinfo <- predictGPdevdevdevdev(fitted_info,t(as.matrix(centercoords[k,])), inputtrain, ftrain)
  predmeandevdevdevdev[k] <- predinfo$pred_mean
  predLBdevdevdevdev[k] = predmeandevdevdevdev[k] - qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
  predUBdevdevdevdev[k] = predmeandevdevdevdev[k]+ qnorm(1-0.05/2)*sqrt(predinfo$pred_var)
}
print(predmeandevdevdevdev[inds])
[1] -0.75730602  3.59103038  0.01582269  1.70411612 -0.96654698 -0.40725528  0.17917950
print(predLBdevdevdevdev[inds])
[1] -4.37970356  3.59102976  0.01582207  1.70411550 -0.96654760 -0.40725590  0.17917888
print(predUBdevdevdevdev[inds])
[1]  2.86509152  3.59103100  0.01582331  1.70411674 -0.96654636 -0.40725466  0.17918012
plot(centercoords[,1],predmeandevdevdevdev, ylim = c(min(predLBdevdevdevdev),max(predUBdevdevdevdev)))
lines(centercoords[,1],predLBdevdevdevdev)
lines(centercoords[,1],predUBdevdevdevdev)

kBT = (293*0.001985)
print((exp(predmean[inds[4]]-min(predmean))/kBT)*sqrt(-2*pi*kBT/predmeandevdev[inds[4]]) + kBT*((1/8)*(predmeandevdevdevdev[inds[4]]/(predmeandevdev[inds[4]])^2) - (5/24)*((predmeandevdevdev[inds[4]])^2/(predmeandevdev[inds[4]])^3)))
[1] 12.20189
print((exp(predLB[inds[4]]-min(predLB))/kBT)*sqrt(-2*pi*kBT/predLBdevdev[inds[4]])+ kBT*((1/8)*(predLBdevdevdevdev[inds[4]]/(predLBdevdev[inds[4]])^2) - (5/24)*((predLBdevdevdev[inds[4]])^2/(predLBdevdev[inds[4]])^3)))
[1] 4.171834
LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQoKYGBge3J9CmxpYnJhcnkocmVhZHIpCmxpYnJhcnkodGFSaWZ4KQpsaWJyYXJ5KHRpZHlyKQpsaWJyYXJ5KHBsb3RseSkKYGBgCmBgYHtyfQpGcmVlRV9Qcm9wYW5lZGlvbC5kZiA8LSByZWFkX3RhYmxlMigiL2hvbWUvYWFyY2hlci9Eb3dubG9hZHMvRnJlZUVfUHJvcGFuZWRpb2wudHh0IikKRnJlZUVfUHJvcGFuYWwuZGYgPC0gcmVhZF90YWJsZTIoIi9ob21lL2FhcmNoZXIvRG93bmxvYWRzL0ZyZWVFX1Byb3BhbmFsWzJdLnR4dCIpCgpGcmVlRV9Qcm9wYW5lZGlvbC5kZiA8LSBGcmVlRV9Qcm9wYW5lZGlvbC5kZltGcmVlRV9Qcm9wYW5lZGlvbC5kZiRgQXhpYWwtY29vcmRgID4gLTkmIEZyZWVFX1Byb3BhbmVkaW9sLmRmJGBBeGlhbC1jb29yZGAgPCA5LF0KRnJlZUVfUHJvcGFuYWwuZGYgPC0gRnJlZUVfUHJvcGFuYWwuZGZbRnJlZUVfUHJvcGFuYWwuZGYkYEF4aWFsLWNvb3JkYCA+IC05JiBGcmVlRV9Qcm9wYW5hbC5kZiRgQXhpYWwtY29vcmRgIDwgOSxdCgpGcmVlRV9Qcm9wYW5lZGlvbF9DZW50ZXIuZGYgPC0gRnJlZUVfUHJvcGFuZWRpb2wuZGZbRnJlZUVfUHJvcGFuZWRpb2wuZGYkYFJhZGlhbC1jb29yZGAgPT0gMC43NSAmIEZyZWVFX1Byb3BhbmVkaW9sLmRmJGBBeGlhbC1jb29yZGAgPiAtOSYgRnJlZUVfUHJvcGFuZWRpb2wuZGYkYEF4aWFsLWNvb3JkYCA8IDksXQpGcmVlRV9Qcm9wYW5hbF9DZW50ZXIuZGYgPC0gRnJlZUVfUHJvcGFuYWwuZGZbRnJlZUVfUHJvcGFuYWwuZGYkYFJhZGlhbC1jb29yZGAgPT0gMC43NSAmIEZyZWVFX1Byb3BhbmFsLmRmJGBBeGlhbC1jb29yZGAgPiAtOSYgRnJlZUVfUHJvcGFuYWwuZGYkYEF4aWFsLWNvb3JkYCA8IDksXQpgYGAKCmBgYHtyfQpnZ3Bsb3QoZGF0YSA9IEZyZWVFX1Byb3BhbmVkaW9sLmRmLCBtYXBwaW5nID0gYWVzKHggPSBgUmFkaWFsLWNvb3JkYCwgeSA9IGBBeGlhbC1jb29yZGAsIGZpbGwgPSBmcmVlRSkpICsgZ2VvbV90aWxlKCkgKyAgeGxhYihsYWJlbCA9ICJSYWRpdXMgQ29vcmQiKSsgIHlsYWIobGFiZWwgPSAiQXhpYWwgQ29vcmQiKSArIGdndGl0bGUoJ0ZyZWUgRW5lcmd5IERpc3RyaWJ1dGlvbiBvZiAxLDItUERPJykKcGxvdChGcmVlRV9Qcm9wYW5lZGlvbF9DZW50ZXIuZGYkYEF4aWFsLWNvb3JkYCxGcmVlRV9Qcm9wYW5lZGlvbF9DZW50ZXIuZGYkZnJlZUUseWxhYj0nQXZlcmFnZSBGcmVlIEVuZXJneSBvZiAxLDItUERPJywgeGxhYj0nQXhpYWwgQ29vcmQnKQpgYGAKYGBge3J9CmdncGxvdChkYXRhID0gRnJlZUVfUHJvcGFuYWwuZGYsIG1hcHBpbmcgPSBhZXMoeCA9IGBSYWRpYWwtY29vcmRgLCB5ID0gYEF4aWFsLWNvb3JkYCwgZmlsbCA9IGZyZWVFKSkgKyBnZW9tX3RpbGUoKSArICB4bGFiKGxhYmVsID0gIlJhZGl1cyBDb29yZCIpKyAgeWxhYihsYWJlbCA9ICJBeGlhbCBDb29yZCIpICsgZ2d0aXRsZSgnRnJlZSBFbmVyZ3kgRGlzdHJpYnV0aW9uIG9mIDEsMi1QRE8nKQpwbG90KEZyZWVFX1Byb3BhbmFsX0NlbnRlci5kZiRgQXhpYWwtY29vcmRgLEZyZWVFX1Byb3BhbmFsX0NlbnRlci5kZiRmcmVlRSx5bGFiPSdBdmVyYWdlIEZyZWUgRW5lcmd5IG9mIFByb3BvbmFsJywgeGxhYj0nQXhpYWwgQ29vcmQnKQpgYGAKCgpgYGB7cn0KY29yciA8LSBmdW5jdGlvbih0djEsdHYyLGV0YTEpewogIEsgPSBtYXRyaXgoMSxucm93PWRpbSh0djEpWzFdLG5jb2w9ZGltKHR2MilbMV0pCiAgZm9yKGRsY3YgaW4gMTpkaW0odHYxKVsyXSl7CiAgSyA8LSBLKmV4cCgtZXRhMVtkbGN2XSpvdXRlcih0djFbLGRsY3ZdLHR2MlssZGxjdl0sJy0nKV4yKQogIH0KICBLCn0KCgoKTUxFIDwtZnVuY3Rpb24obG9nZXRhLCB0bywgZm8pewogIG0gPC0gbGVuZ3RoKGxvZ2V0YSkKICBsb2dldGExIDwtIGxvZ2V0YVstbV0KICBsb2dldGEyIDwtIGxvZ2V0YVttXQogIAogICMgb2J0YWluIHRoZSBjb3ZhcmlhbmNlIHZpYSBrcm9uZWNrZXIgcHJvZHVjdAogIAogIGNvcnIwIDwtIGNvcnIodG8sIHRvLCBleHAobG9nZXRhMSkpIyBjKC4sIC4pCiAga2FwcGEwIDwtIGNvcnIwICsgZXhwKGxvZ2V0YTIpKmRpYWcobGVuZ3RoKGFzLnZlY3RvcihmbykpKQogIGthcHBhMF9pbnYgPC0gc29sdmUoa2FwcGEwKQogIGthcHBhMF9pbnZfZjAgPC0gc29sdmUoa2FwcGEwLCBhcy52ZWN0b3IoZm8pKQogIGthcHBhMF9pbnZfMSA8LSBzb2x2ZShrYXBwYTAsIHJlcCgxLCBsZW5ndGgoYXMudmVjdG9yKGZvKSkpKSMgZXN0aW1hdGVkIHZhbHVlcwogIGdhbW1haGF0IDwtIHN1bShrYXBwYTBfaW52X2YwKS9zdW0oa2FwcGEwX2ludl8xKQogIHNpZ21haGF0IDwtIG1lYW4odChrYXBwYTBfaW52X2YwIC0gZ2FtbWFoYXQqa2FwcGEwX2ludl8xKSooYXMudmVjdG9yKGZvKSAtIGdhbW1haGF0KSkKICByZXR1cm4obGlzdChnYW1tYWhhdD1nYW1tYWhhdCwgc2lnbWFoYXQ9c2lnbWFoYXQpKQogIH0KYGBgCgpgYGB7cn0KIyB0cmFpbiBlYWNoIGYgYXMgR1AKbmVnbG9nbGlrIDwtZnVuY3Rpb24obG9nZXRhLCB0MCwgZjApewogIG4gPC0gbGVuZ3RoKGFzLnZlY3RvcihmMCkpIyBvYnRhaW4gc2lnbWFoYXQgZm9yIGEgZ2l2ZW4gZXRhCiAgbSA8LSBsZW5ndGgobG9nZXRhKQogIGxvZ2V0YTEgPC0gbG9nZXRhWy1tXQogIGxvZ2V0YTIgPC0gbG9nZXRhW21dCiAgZXN0aW1hdG9ycyA8LSBNTEUobG9nZXRhLCB0MCwgZjApCiAgc2lnbWFoYXQgPC0gZXN0aW1hdG9ycyRzaWdtYWhhdCMgb2J0YWluIHRoZSBjb3ZhcmlhbmNlIG1hdHJpY2VzCiAgY29ycjAgPC0gY29ycih0MCwgdDAsIGV4cChsb2dldGExKSkjIGModDAsIHQwKQogIGthcHBhMCA8LSBjb3JyMCArIGV4cChsb2dldGEyKSpkaWFnKG4pCgogICMgb2J0YWluIHRoZSBsb2cgb2YgZGV0ZXJtaW5hbnRzIGVmZmljaWVudGx5CiAgbG9nZGV0a2FwcGEwIDwtIGRldGVybWluYW50KGthcHBhMCxsb2dhcml0aG09VFJVRSkkbW9kdWx1cwogIHJldHVybigxLzIqbG9nZGV0a2FwcGEwICsgbi8yKmxvZyhzaWdtYWhhdCkpCiAgfQpgYGAKCmBgYHtyfQpmaXRHUCA8LWZ1bmN0aW9uKGxvZ2V0YT1jKDEsIDEsIDEpLCBsb3dlcmI9YygtMSwgLTEsLTEpLCB1cHBlcmI9Yyg0LCA0LDQpLCB0X3RyLCBmX3RyKXsKICBsb2dfZXRhX2hhdCA8LSBvcHRpbShwYXI9bG9nZXRhLGZuPW5lZ2xvZ2xpayxsb3dlcj1sb3dlcmIsdXBwZXI9dXBwZXJiLG1ldGhvZD0iTC1CRkdTLUIiLCB0MD10X3RyLCBmMD1mX3RyKSRwYXIKICBlc3RpbWF0b3JzIDwtIE1MRShsb2dfZXRhX2hhdCwgdF90ciwgZl90cikKICBtIDwtIGxlbmd0aChsb2dfZXRhX2hhdCkKICByZXR1cm4obGlzdChldGFoYXQ9ZXhwKGxvZ19ldGFfaGF0Wy1tXSksZz1leHAobG9nX2V0YV9oYXRbbV0pLHNpZ21haGF0PWVzdGltYXRvcnMkc2lnbWFoYXQsZ2FtbWFoYXQ9ZXN0aW1hdG9ycyRnYW1tYWhhdCkpCn0KcHJlZGljdEdQIDwtZnVuY3Rpb24oZml0dGVkX2luZm8sIHRfdGVzdCwgdF90ciwgZl90cil7IyBvYnRhaW4gdGhlIGZpdHRlZCBpbmZvCiAgZXRhaGF0MSA8LSBmaXR0ZWRfaW5mbyRldGFoYXQKICBnIDwtIGZpdHRlZF9pbmZvJGcKICBzaWdtYWhhdCA8LSBmaXR0ZWRfaW5mbyRzaWdtYWhhdAogIGdhbW1haGF0IDwtIGZpdHRlZF9pbmZvJGdhbW1haGF0CiAgCiAgIyBvYnRhaW4gdGhlIGZpbmFsIGNvdmFyaWFuY2Ugc3RydWN0dXJlCiAgCiAgY29ycjAgPC0gY29ycih0X3RyLCB0X3RyLCBldGFoYXQxKSAgIyBjKHRfdHIsIHRfdHIpCiAga2FwcGEwIDwtIGNvcnIwICsgZypkaWFnKGxlbmd0aChhcy52ZWN0b3IoZl90cikpKQogIGthcHBhMF9pbnYgPC0gKDEvc2lnbWFoYXQpKnNvbHZlKGthcHBhMCkKICBrYXBwYTBfaW52X2YwIDwtIGthcHBhMF9pbnYlKiUoYXMudmVjdG9yKGZfdHIpIC0gZ2FtbWFoYXQpCiAgY29ycl9UdiA8LSBjb3JyKHRfdGVzdCwgdF90ciwgZXRhaGF0MSkjIGModF90ZXN0LCB0X3RyKQogIGthcHBhX1R2IDwtIHNpZ21haGF0KmNvcnJfVHYjIGModF90ZXN0LCB0X3RyKQogIAogIGNvcnJfdG9UdiA8LSBjb3JyKHRfdHIsIHRfdGVzdCwgZXRhaGF0MSkjIGModF90ciwgdF90ZXN0KQogIGthcHBhX3RvVHYgPC0gc2lnbWFoYXQgKmNvcnJfdG9UdiMgIGMoLiwgLikKICAKICBnYW1tYV92ZWMgPC0gKGdhbW1haGF0ICsga2FwcGFfVHYlKiVrYXBwYTBfaW52X2YwKQogIGRrYXBwYSA8LSBzaWdtYWhhdCAqIHJlcCgxLCBkaW0odF90ZXN0KVsxXSpkaW0oZl90cilbMl0pCiAga2FwcGFfdmVjIDwtIChwbWF4KGRrYXBwYSAtYXBwbHkodCgoa2FwcGFfVHYpJSolKGthcHBhMF9pbnYpKSooa2FwcGFfdG9UdiksIDIsIHN1bSksIDEwKiooLTEzKSkpCiAgcmV0dXJuKGxpc3QocHJlZF9tZWFuPWdhbW1hX3ZlYywgcHJlZF92YXI9a2FwcGFfdmVjKSl9CmBgYAoKYGBge3J9CgpmdHJhaW4gPC0gYXMubWF0cml4KEZyZWVFX1Byb3BhbmFsLmRmWywzXSkKaW5wdXR0cmFpbiA8LWFzLm1hdHJpeChGcmVlRV9Qcm9wYW5hbC5kZlssMToyXSkKZml0dGVkX2luZm8gPC0gZml0R1AobG9nZXRhPWMoMSwgMSwxKSwgbG93ZXJiPWMoMCwgMCwwKSwgdXBwZXJiPWMoNSwgNSw1KSwgaW5wdXR0cmFpbiwgZnRyYWluKQpmaXR0ZWRfaW5mbwpgYGAKCgpgYGB7cn0KcHJlZG1lYW4gPXJlcCgwLGRpbShpbnB1dHRyYWluKVsxXSkKcHJlZExCID0gIHJlcCgwLGRpbShpbnB1dHRyYWluKVsxXSkKcHJlZFVCID1yZXAoMCxkaW0oaW5wdXR0cmFpbilbMV0pCmZvcihrIGluIDE6ZGltKGlucHV0dHJhaW4pWzFdKXsKICBwcmVkaW5mbyA8LSBwcmVkaWN0R1AoZml0dGVkX2luZm8sdChhcy5tYXRyaXgoaW5wdXR0cmFpbltrLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuW2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQltrXSA9IHByZWRtZWFuW2tdLSBxbm9ybSgxLTAuMDUvMikqc3FydChwcmVkaW5mbyRwcmVkX3ZhcikKICBwcmVkVUJba10gPSBwcmVkbWVhbltrXSsgcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCn0KcHJpbnQocGFzdGUoJ1BDR1BNZXRob2Q6IFRlc3QgTVNFOicscm91bmQobWVhbigocHJlZG1lYW4tZnRyYWluKV4yKSwzKSwnLCBDb3ZlcmFnZTonLDEwMCpyb3VuZChtZWFuKChwcmVkTEI8ZnRyYWluKSoocHJlZFVCPmZ0cmFpbikpLDMpLCclJyxzZXAgPScnKSkKYGBgCmBgYHtyfQpGcmVlRV9Qcm9wYW5hbF9QcmVkLmRmIDwtIEZyZWVFX1Byb3BhbmFsLmRmCkZyZWVFX1Byb3BhbmFsX0xCLmRmIDwtIEZyZWVFX1Byb3BhbmFsLmRmCkZyZWVFX1Byb3BhbmFsX1VCLmRmIDwtIEZyZWVFX1Byb3BhbmFsLmRmCgoKRnJlZUVfUHJvcGFuYWxfUHJlZC5kZiRmcmVlRSA8LSBwcmVkbWVhbgpGcmVlRV9Qcm9wYW5hbF9MQi5kZiRmcmVlRSA8LSBwcmVkTEIKRnJlZUVfUHJvcGFuYWxfVUIuZGYkZnJlZUUgPC0gcHJlZFVCCnogPC0gYXMubWF0cml4KHNwcmVhZChGcmVlRV9Qcm9wYW5hbC5kZiwgYFJhZGlhbC1jb29yZGAsIGZyZWVFKVssLTFdKQp6cHJlZCA8LSBhcy5tYXRyaXgoc3ByZWFkKEZyZWVFX1Byb3BhbmFsX1ByZWQuZGYsIGBSYWRpYWwtY29vcmRgLCBmcmVlRSlbLC0xXSkKekxCIDwtIGFzLm1hdHJpeChzcHJlYWQoRnJlZUVfUHJvcGFuYWxfTEIuZGYsIGBSYWRpYWwtY29vcmRgLCBmcmVlRSlbLC0xXSkKelVCIDwtIGFzLm1hdHJpeChzcHJlYWQoRnJlZUVfUHJvcGFuYWxfVUIuZGYsIGBSYWRpYWwtY29vcmRgLCBmcmVlRSlbLC0xXSkKCmNvbG9yIDwtIHJlcCgxLCBsZW5ndGgoeikpCmRpbShjb2xvcikgPC0gZGltKHopCgp4IDwtdW5pcXVlKEZyZWVFX1Byb3BhbmFsLmRmJGBSYWRpYWwtY29vcmRgKQp5IDwtdW5pcXVlKEZyZWVFX1Byb3BhbmFsLmRmJGBBeGlhbC1jb29yZGApCgpwPC0gIHBsb3RfbHkoKSAlPiUgYWRkX3RyYWNlKHg9RnJlZUVfUHJvcGFuYWwuZGYkYFJhZGlhbC1jb29yZGAsIHk9RnJlZUVfUHJvcGFuYWwuZGYkYEF4aWFsLWNvb3JkYCwgej1GcmVlRV9Qcm9wYW5hbC5kZiRmcmVlRSwgbW9kZT0gIm1hcmtlcnMiLHR5cGUgPSAic2NhdHRlcjNkIixtYXJrZXIgPSBsaXN0KHNpemUgPSAyLCBjb2xvciA9ICJibHVlIiwgc3ltYm9sID0gMTA0KSwgb3BhY2l0eT0wLjcsc2hvd2xlZ2VuZCA9IEZBTFNFICkKcCA8LSBwICU+JSBhZGRfdHJhY2UoIHg9eCwgeT15LCB6PXpwcmVkLHR5cGUgPSAic3VyZmFjZSIsY29sb3JzY2FsZSA9IGxpc3QoYygwLCAxKSwgYygidGFuIiwgImJsdWUiKSksb3BhY2l0eT0wLjcpICAKCmRmMSA8LSBzcGxpdChGcmVlRV9Qcm9wYW5hbF9MQi5kZiwgRnJlZUVfUHJvcGFuYWxfTEIuZGYkYFJhZGlhbC1jb29yZGApCmRmMiA8LSBzcGxpdChGcmVlRV9Qcm9wYW5hbF9MQi5kZiwgRnJlZUVfUHJvcGFuYWxfTEIuZGYkYEF4aWFsLWNvb3JkYCkKCmZvcihpIGluIHNlcV9hbG9uZyhkZjEpKXsKICBkZl9zcCA8LSBkZjFbW2ldXQogIHAgPC0gYWRkX3RyYWNlKHAsIGxpbmUgPSBsaXN0KAogICAgY29sb3IgPSAiIzAwMDAwMCIsIAogICAgd2lkdGggPSAyCiAgKSwgCiAgbW9kZSA9ICJsaW5lcyIsIAogIHR5cGUgPSAic2NhdHRlcjNkIiwgCiAgeCA9IGRmX3NwJGBSYWRpYWwtY29vcmRgLAogIHkgPSBkZl9zcCRgQXhpYWwtY29vcmRgLAogIHogPSBkZl9zcCRmcmVlRSwKICBzaG93bGVnZW5kID0gRkFMU0UpCn0KZm9yKGkgaW4gc2VxX2Fsb25nKGRmMikpewogIGRmX3NwIDwtIGRmMltbaV1dCiAgcCA8LSBhZGRfdHJhY2UocCwgbGluZSA9IGxpc3QoCiAgICBjb2xvciA9ICIjMDAwMDAwIiwgCiAgICB3aWR0aCA9IDIKICApLCAKICBtb2RlID0gImxpbmVzIiwgCiAgdHlwZSA9ICJzY2F0dGVyM2QiLCAKICB4ID0gZGZfc3AkYFJhZGlhbC1jb29yZGAsCiAgeSA9IGRmX3NwJGBBeGlhbC1jb29yZGAsCiAgeiA9IGRmX3NwJGZyZWVFLAogIHNob3dsZWdlbmQgPSBGQUxTRSkKfQoKZGYxIDwtIHNwbGl0KEZyZWVFX1Byb3BhbmFsX1VCLmRmLCBGcmVlRV9Qcm9wYW5hbF9VQi5kZiRgUmFkaWFsLWNvb3JkYCkKZGYyIDwtIHNwbGl0KEZyZWVFX1Byb3BhbmFsX1VCLmRmLCBGcmVlRV9Qcm9wYW5hbF9VQi5kZiRgQXhpYWwtY29vcmRgKQoKZm9yKGkgaW4gc2VxX2Fsb25nKGRmMSkpewogIGRmX3NwIDwtIGRmMVtbaV1dCiAgcCA8LSBhZGRfdHJhY2UocCwgbGluZSA9IGxpc3QoCiAgICBjb2xvciA9ICIjRkYwMDAwIiwgCiAgICB3aWR0aCA9IDIKICApLCAKICBtb2RlID0gImxpbmVzIiwgCiAgdHlwZSA9ICJzY2F0dGVyM2QiLCAKICB4ID0gZGZfc3AkYFJhZGlhbC1jb29yZGAsCiAgeSA9IGRmX3NwJGBBeGlhbC1jb29yZGAsCiAgeiA9IGRmX3NwJGZyZWVFLAogIHNob3dsZWdlbmQgPSBGQUxTRSkKfQpmb3IoaSBpbiBzZXFfYWxvbmcoZGYyKSl7CiAgZGZfc3AgPC0gZGYyW1tpXV0KICBwIDwtIGFkZF90cmFjZShwLCBsaW5lID0gbGlzdCgKICAgIGNvbG9yID0gIiNGRjAwMDAiLCAKICAgIHdpZHRoID0gMgogICksIAogIG1vZGUgPSAibGluZXMiLCAKICB0eXBlID0gInNjYXR0ZXIzZCIsIAogIHggPSBkZl9zcCRgUmFkaWFsLWNvb3JkYCwKICB5ID0gZGZfc3AkYEF4aWFsLWNvb3JkYCwKICB6ID0gZGZfc3AkZnJlZUUsCiAgc2hvd2xlZ2VuZCA9IEZBTFNFKQp9CnAlPiUgbGF5b3V0KHNjZW5lID0gbGlzdCggeGF4aXMgPSBsaXN0KHRpdGxlID0gIlJhZGl1cywgciIpLHlheGlzID0gbGlzdCh0aXRsZSA9ICJBeGlhbCwgeiIpKSkKYGBgCgpgYGB7cn0KCmF4aWFsbGVuIDwtIDEwMDAKYXhpYWxjb29yZCA8LSB1bmlxdWUoRnJlZUVfUHJvcGFuYWwuZGYkYEF4aWFsLWNvb3JkYCkKY2VudGVyY29vcmRzIDwtYXMubWF0cml4KCBkYXRhLmZyYW1lKHNlcShtaW4oYXhpYWxjb29yZCksbWF4KGF4aWFsY29vcmQpLGxlbmd0aC5vdXQ9YXhpYWxsZW4pLHJlcCgwLGF4aWFsbGVuKSkpCgpwcmVkbWVhbiA9cmVwKDAsZGltKGNlbnRlcmNvb3JkcylbMV0pCnByZWRMQiA9ICByZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZFVCID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKZm9yKGsgaW4gMTpkaW0oY2VudGVyY29vcmRzKVsxXSl7CiAgcHJlZGluZm8gPC0gcHJlZGljdEdQKGZpdHRlZF9pbmZvLHQoYXMubWF0cml4KGNlbnRlcmNvb3Jkc1trLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuW2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQltrXSA9IHByZWRtZWFuW2tdLSBxbm9ybSgxLTAuMDUvMikqc3FydChwcmVkaW5mbyRwcmVkX3ZhcikKICBwcmVkVUJba10gPSBwcmVkbWVhbltrXSsgcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCn0KYGBgCgpgYGB7cn0KcGxvdChjZW50ZXJjb29yZHNbLDFdLHByZWRtZWFuLCB5bGltID0gYyhtaW4ocHJlZExCKSxtYXgocHJlZFVCKSkpCmxpbmVzKGNlbnRlcmNvb3Jkc1ssMV0scHJlZExCKQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRVQikKZGVsdGF4ID0gKG1heChheGlhbGNvb3JkKS1taW4oYXhpYWxjb29yZCkpLyhheGlhbGxlbi0xKQpwcmVkbWVhbmRldiA9IChwcmVkbWVhblszOmF4aWFsbGVuXS1wcmVkbWVhblsxOihheGlhbGxlbi0yKV0pL2RlbHRheApwcmVkZGV2TEIgPSAocHJlZExCWzM6YXhpYWxsZW5dLXByZWRMQlsxOihheGlhbGxlbi0yKV0pL2RlbHRheApwcmVkZGV2VUIgPSAocHJlZFVCWzM6YXhpYWxsZW5dLXByZWRVQlsxOihheGlhbGxlbi0yKV0pL2RlbHRheApwcmVkbWVhbmRldmRldiA9IChwcmVkbWVhblszOmF4aWFsbGVuXS0yKnByZWRtZWFuWzI6KGF4aWFsbGVuLTEpXSArIHByZWRtZWFuWzE6KGF4aWFsbGVuLTIpXSkvZGVsdGF4XjIKcHJlZGRldmRldkxCID0gKHByZWRMQlszOmF4aWFsbGVuXS0yKnByZWRtZWFuWzI6KGF4aWFsbGVuLTEpXSArcHJlZExCWzE6KGF4aWFsbGVuLTIpXSkvZGVsdGF4XjIKcHJlZGRldmRldlVCID0gKHByZWRVQlszOmF4aWFsbGVuXS0yKnByZWRtZWFuWzI6KGF4aWFsbGVuLTEpXSArcHJlZFVCWzE6KGF4aWFsbGVuLTIpXSkvZGVsdGF4XjIKCmluZG1lYW4gPSB3aGljaChwcmVkbWVhbj09bWF4KHByZWRtZWFuKSkKaW5kTEIgPSB3aGljaChwcmVkTEI9PW1heChwcmVkTEIpKQppbmRVQiA9IHdoaWNoKHByZWRVQj09bWF4KHByZWRVQikpCgpwcmludChwcmVkbWVhbltpbmRtZWFuXSkKcHJpbnQocHJlZExCW2luZExCXSkKcHJpbnQocHJlZFVCW2luZFVCXSkKCnBsb3QoY2VudGVyY29vcmRzWzI6KGF4aWFsbGVuLTEpLDFdLHByZWRtZWFuZGV2ICwgeWxpbSA9IGMobWluKHByZWRtZWFuZGV2LHByZWRkZXZMQixwcmVkZGV2VUIpLG1heChwcmVkbWVhbmRldixwcmVkZGV2TEIscHJlZGRldlVCKSkpCmxpbmVzKGNlbnRlcmNvb3Jkc1syOihheGlhbGxlbi0xKSwxXSxwcmVkZGV2TEIpCmxpbmVzKGNlbnRlcmNvb3Jkc1syOihheGlhbGxlbi0xKSwxXSxwcmVkZGV2VUIpCgpwcmludChwcmVkbWVhbmRldltpbmRtZWFuXSkKcHJpbnQocHJlZGRldkxCW2luZExCXSkKcHJpbnQocHJlZGRldlVCW2luZFVCXSkKCnBsb3QoY2VudGVyY29vcmRzWzI6KGF4aWFsbGVuLTEpLDFdLHByZWRtZWFuZGV2ZGV2ICwgeWxpbSA9IGMobWluKHByZWRtZWFuZGV2ZGV2LHByZWRkZXZkZXZMQixwcmVkZGV2ZGV2VUIpLG1heChwcmVkbWVhbmRldmRldixwcmVkZGV2ZGV2TEIscHJlZGRldmRldlVCKSkpCmxpbmVzKGNlbnRlcmNvb3Jkc1syOihheGlhbGxlbi0xKSwxXSxwcmVkZGV2ZGV2TEIpCmxpbmVzKGNlbnRlcmNvb3Jkc1syOihheGlhbGxlbi0xKSwxXSxwcmVkZGV2ZGV2VUIpCgpwcmludChwcmVkbWVhbmRldmRldltpbmRtZWFuXSkKcHJpbnQocHJlZGRldmRldkxCW2luZExCXSkKcHJpbnQocHJlZGRldmRldlVCW2luZFVCXSkKCgpgYGAKCmBgYHtyfQpjb3JyZGV2b2J2IDwtIGZ1bmN0aW9uKHR2MSx0djIsZXRhMSl7CiAgSyA8LSAtMipldGExWzFdKmV4cCgtZXRhMVsxXSpvdXRlcih0djFbLDFdLHR2MlssMV0sJy0nKV4yKSoob3V0ZXIodHYxWywxXSx0djJbLDFdLCctJykpCiAgSyA8LSBLKmV4cCgtZXRhMVsyXSphYnMob3V0ZXIodHYxWywyXSx0djJbLDJdLCctJykpXjIpCiAgSwp9CgpwcmVkaWN0R1BkZXYgPC1mdW5jdGlvbihmaXR0ZWRfaW5mbywgdF90ZXN0LCB0X3RyLCBmX3RyKXsjIG9idGFpbiB0aGUgZml0dGVkIGluZm8KICBldGFoYXQxIDwtIGZpdHRlZF9pbmZvJGV0YWhhdAogIGcgPC0gZml0dGVkX2luZm8kZwogIHNpZ21haGF0IDwtIGZpdHRlZF9pbmZvJHNpZ21haGF0CiAgZ2FtbWFoYXQgPC0gZml0dGVkX2luZm8kZ2FtbWFoYXQKICAKICAjIG9idGFpbiB0aGUgZmluYWwgY292YXJpYW5jZSBzdHJ1Y3R1cmUKICAKICBjb3JyMCA8LSBjb3JyKHRfdHIsIHRfdHIsIGV0YWhhdDEpICAjIGModF90ciwgdF90cikKICBrYXBwYTAgPC0gY29ycjAgKyBnKmRpYWcobGVuZ3RoKGFzLnZlY3RvcihmX3RyKSkpCiAga2FwcGEwX2ludiA8LSAoMS9zaWdtYWhhdCkqc29sdmUoa2FwcGEwKQogIGthcHBhMF9pbnZfZjAgPC0ga2FwcGEwX2ludiUqJShhcy52ZWN0b3IoZl90cikgLSBnYW1tYWhhdCkKICBjb3JyX1R2IDwtIGNvcnJkZXZvYnYodF90ZXN0LCB0X3RyLCBldGFoYXQxKSMgYyh0X3Rlc3QsIHRfdHIpCiAga2FwcGFfVHYgPC0gc2lnbWFoYXQqY29ycl9UdiMgYyh0X3Rlc3QsIHRfdHIpCiAgCiAgY29ycl90b1R2IDwtIGNvcnJkZXZvYnYodF90ciwgdF90ZXN0LCBldGFoYXQxKSMgYyh0X3RyLCB0X3Rlc3QpCiAga2FwcGFfdG9UdiA8LSBzaWdtYWhhdCAqY29ycl90b1R2IyAgYyguLCAuKQogIAogIGdhbW1hX3ZlYyA8LSAga2FwcGFfVHYlKiVrYXBwYTBfaW52X2YwCiAgZGthcHBhIDwtIHNpZ21haGF0KigxICsgZykqcmVwKDEsIGRpbSh0X3Rlc3QpWzFdKmRpbShmX3RyKVsyXSkKICBrYXBwYV92ZWMgPC0gKHBtYXgoZGthcHBhIC1hcHBseSh0KChrYXBwYV9UdiklKiUoa2FwcGEwX2ludikpKihrYXBwYV90b1R2KSwgMiwgc3VtKSwgMTAqKigtMTMpKSkKICByZXR1cm4obGlzdChwcmVkX21lYW49Z2FtbWFfdmVjLCBwcmVkX3Zhcj1rYXBwYV92ZWMpKX0KYGBgCgpgYGB7cn0KCnByZWRtZWFuZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZExCZGV2ID0gIHJlcCgwLGRpbShjZW50ZXJjb29yZHMpWzFdKQpwcmVkVUJkZXYgPXJlcCgwLGRpbShjZW50ZXJjb29yZHMpWzFdKQpmb3IoayBpbiAxOmRpbShjZW50ZXJjb29yZHMpWzFdKXsKICBwcmVkaW5mbyA8LSBwcmVkaWN0R1BkZXYoZml0dGVkX2luZm8sdChhcy5tYXRyaXgoY2VudGVyY29vcmRzW2ssXSkpLCBpbnB1dHRyYWluLCBmdHJhaW4pCiAgcHJlZG1lYW5kZXZba10gPC0gcHJlZGluZm8kcHJlZF9tZWFuCiAgcHJlZExCZGV2W2tdID0gcHJlZG1lYW5kZXZba10gLSBxbm9ybSgxLTAuMDUvMikqc3FydChwcmVkaW5mbyRwcmVkX3ZhcikKICBwcmVkVUJkZXZba10gPSBwcmVkbWVhbmRldltrXSsgcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCn0KYGBgCgoKYGBge3J9CnBsb3QoY2VudGVyY29vcmRzWywxXSxwcmVkbWVhbmRldikjIHlsaW0gPSBjKG1pbihwcmVkTEJkZXYpLG1heChwcmVkVUJkZXYpKSkKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkTEJkZXYpCmxpbmVzKGNlbnRlcmNvb3Jkc1ssMV0scHJlZFVCZGV2KQppbmRzID0gd2hpY2goZGlmZihzaWduKHByZWRtZWFuZGV2KSkgIT0gMCkKcHJpbnQocHJlZG1lYW5kZXZbaW5kc10pCnByaW50KHByZWRMQmRldltpbmRzXSkKcHJpbnQocHJlZFVCZGV2W2luZHNdKQpgYGAKCmBgYHtyfQpjb3JyZGV2ZGV2b2J2IDwtIGZ1bmN0aW9uKHR2MSx0djIsZXRhMSl7CiAgSyA8LSAyKmV0YTFbMV0qZXhwKC1ldGExWzFdKm91dGVyKHR2MVssMV0sdHYyWywxXSwnLScpXjIpKigyKmV0YTFbMV0qKG91dGVyKHR2MVssMV0sdHYyWywxXSwnLScpKV4yLTEpCiAgSyA8LSBLKmV4cCgtZXRhMVsyXSphYnMob3V0ZXIodHYxWywyXSx0djJbLDJdLCctJykpXjIpCiAgSwp9Cgpjb3JyZGV2ZGV2ZGV2b2J2IDwtIGZ1bmN0aW9uKHR2MSx0djIsZXRhMSl7CiAgSyA8LSAtNCooZXRhMVsxXV4yKSoob3V0ZXIodHYxWywxXSx0djJbLDFdLCctJykpKmV4cCgtZXRhMVsxXSpvdXRlcih0djFbLDFdLHR2MlssMV0sJy0nKV4yKSooMipldGExWzFdKihvdXRlcih0djFbLDFdLHR2MlssMV0sJy0nKSleMi0zKQogIEsgPC0gSypleHAoLWV0YTFbMl0qYWJzKG91dGVyKHR2MVssMl0sdHYyWywyXSwnLScpKV4yKQogIEsKfQoKY29ycmRldmRldmRldmRldm9idiA8LSBmdW5jdGlvbih0djEsdHYyLGV0YTEpewogIEsgPC0gNCooZXRhMVsxXV4yKSpleHAoLWV0YTFbMV0qb3V0ZXIodHYxWywxXSx0djJbLDFdLCctJyleMikqKDQqZXRhMVsxXSoob3V0ZXIodHYxWywxXSx0djJbLDFdLCctJyleMikqKGV0YTFbMV0qKG91dGVyKHR2MVssMV0sdHYyWywxXSwnLScpKV4yLTMpKzMpCiAgSyA8LSBLKmV4cCgtZXRhMVsyXSphYnMob3V0ZXIodHYxWywyXSx0djJbLDJdLCctJykpXjIpCiAgSwp9CgpwcmVkaWN0R1BkZXZkZXYgPC1mdW5jdGlvbihmaXR0ZWRfaW5mbywgdF90ZXN0LCB0X3RyLCBmX3RyKXsjIG9idGFpbiB0aGUgZml0dGVkIGluZm8KICBldGFoYXQxIDwtIGZpdHRlZF9pbmZvJGV0YWhhdAogIGcgPC0gZml0dGVkX2luZm8kZwogIHNpZ21haGF0IDwtIGZpdHRlZF9pbmZvJHNpZ21haGF0CiAgZ2FtbWFoYXQgPC0gZml0dGVkX2luZm8kZ2FtbWFoYXQKICAKICAjIG9idGFpbiB0aGUgZmluYWwgY292YXJpYW5jZSBzdHJ1Y3R1cmUKICAKICBjb3JyMCA8LSBjb3JyKHRfdHIsIHRfdHIsIGV0YWhhdDEpICAjIGModF90ciwgdF90cikKICBrYXBwYTAgPC0gY29ycjAgKyBnKmRpYWcobGVuZ3RoKGFzLnZlY3RvcihmX3RyKSkpCiAga2FwcGEwX2ludiA8LSAoMS9zaWdtYWhhdCkqc29sdmUoa2FwcGEwKQogIGthcHBhMF9pbnZfZjAgPC0ga2FwcGEwX2ludiUqJShhcy52ZWN0b3IoZl90cikgLSBnYW1tYWhhdCkKICBjb3JyX1R2IDwtIGNvcnJkZXZkZXZvYnYodF90ZXN0LCB0X3RyLCBldGFoYXQxKSMgYyh0X3Rlc3QsIHRfdHIpCiAga2FwcGFfVHYgPC0gc2lnbWFoYXQqY29ycl9UdiMgYyh0X3Rlc3QsIHRfdHIpCiAgCiAgY29ycl90b1R2IDwtIGNvcnJkZXZkZXZvYnYodF90ciwgdF90ZXN0LCBldGFoYXQxKSMgYyh0X3RyLCB0X3Rlc3QpCiAga2FwcGFfdG9UdiA8LSBzaWdtYWhhdCAqY29ycl90b1R2IyAgYyguLCAuKQogIAogIGdhbW1hX3ZlYyA8LSAga2FwcGFfVHYlKiVrYXBwYTBfaW52X2YwCiAgZGthcHBhIDwtIHNpZ21haGF0KigxICsgZykqcmVwKDEsIGRpbSh0X3Rlc3QpWzFdKmRpbShmX3RyKVsyXSkKICBrYXBwYV92ZWMgPC0gKHBtYXgoZGthcHBhIC1hcHBseSh0KChrYXBwYV9UdiklKiUoa2FwcGEwX2ludikpKihrYXBwYV90b1R2KSwgMiwgc3VtKSwgMTAqKigtMTMpKSkKICByZXR1cm4obGlzdChwcmVkX21lYW49Z2FtbWFfdmVjLCBwcmVkX3Zhcj1rYXBwYV92ZWMpKX0KCnByZWRpY3RHUGRldmRldmRldiA8LWZ1bmN0aW9uKGZpdHRlZF9pbmZvLCB0X3Rlc3QsIHRfdHIsIGZfdHIpeyMgb2J0YWluIHRoZSBmaXR0ZWQgaW5mbwogIGV0YWhhdDEgPC0gZml0dGVkX2luZm8kZXRhaGF0CiAgZyA8LSBmaXR0ZWRfaW5mbyRnCiAgc2lnbWFoYXQgPC0gZml0dGVkX2luZm8kc2lnbWFoYXQKICBnYW1tYWhhdCA8LSBmaXR0ZWRfaW5mbyRnYW1tYWhhdAogIAogICMgb2J0YWluIHRoZSBmaW5hbCBjb3ZhcmlhbmNlIHN0cnVjdHVyZQogIAogIGNvcnIwIDwtIGNvcnIodF90ciwgdF90ciwgZXRhaGF0MSkgICMgYyh0X3RyLCB0X3RyKQogIGthcHBhMCA8LSBjb3JyMCArIGcqZGlhZyhsZW5ndGgoYXMudmVjdG9yKGZfdHIpKSkKICBrYXBwYTBfaW52IDwtICgxL3NpZ21haGF0KSpzb2x2ZShrYXBwYTApCiAga2FwcGEwX2ludl9mMCA8LSBrYXBwYTBfaW52JSolKGFzLnZlY3RvcihmX3RyKSAtIGdhbW1haGF0KQogIGNvcnJfVHYgPC0gY29ycmRldmRldmRldm9idih0X3Rlc3QsIHRfdHIsIGV0YWhhdDEpIyBjKHRfdGVzdCwgdF90cikKICBrYXBwYV9UdiA8LSBzaWdtYWhhdCpjb3JyX1R2IyBjKHRfdGVzdCwgdF90cikKICAKICBjb3JyX3RvVHYgPC0gY29ycmRldmRldmRldm9idih0X3RyLCB0X3Rlc3QsIGV0YWhhdDEpIyBjKHRfdHIsIHRfdGVzdCkKICBrYXBwYV90b1R2IDwtIHNpZ21haGF0ICpjb3JyX3RvVHYjICBjKC4sIC4pCiAgCiAgZ2FtbWFfdmVjIDwtICBrYXBwYV9UdiUqJWthcHBhMF9pbnZfZjAKICBka2FwcGEgPC0gc2lnbWFoYXQqKDEgKyBnKSpyZXAoMSwgZGltKHRfdGVzdClbMV0qZGltKGZfdHIpWzJdKQogIGthcHBhX3ZlYyA8LSAocG1heChka2FwcGEgLWFwcGx5KHQoKGthcHBhX1R2KSUqJShrYXBwYTBfaW52KSkqKGthcHBhX3RvVHYpLCAyLCBzdW0pLCAxMCoqKC0xMykpKQogIHJldHVybihsaXN0KHByZWRfbWVhbj1nYW1tYV92ZWMsIHByZWRfdmFyPWthcHBhX3ZlYykpfQoKcHJlZGljdEdQZGV2ZGV2ZGV2ZGV2IDwtZnVuY3Rpb24oZml0dGVkX2luZm8sIHRfdGVzdCwgdF90ciwgZl90cil7IyBvYnRhaW4gdGhlIGZpdHRlZCBpbmZvCiAgZXRhaGF0MSA8LSBmaXR0ZWRfaW5mbyRldGFoYXQKICBnIDwtIGZpdHRlZF9pbmZvJGcKICBzaWdtYWhhdCA8LSBmaXR0ZWRfaW5mbyRzaWdtYWhhdAogIGdhbW1haGF0IDwtIGZpdHRlZF9pbmZvJGdhbW1haGF0CiAgCiAgIyBvYnRhaW4gdGhlIGZpbmFsIGNvdmFyaWFuY2Ugc3RydWN0dXJlCiAgCiAgY29ycjAgPC0gY29ycih0X3RyLCB0X3RyLCBldGFoYXQxKSAgIyBjKHRfdHIsIHRfdHIpCiAga2FwcGEwIDwtIGNvcnIwICsgZypkaWFnKGxlbmd0aChhcy52ZWN0b3IoZl90cikpKQogIGthcHBhMF9pbnYgPC0gKDEvc2lnbWFoYXQpKnNvbHZlKGthcHBhMCkKICBrYXBwYTBfaW52X2YwIDwtIGthcHBhMF9pbnYlKiUoYXMudmVjdG9yKGZfdHIpIC0gZ2FtbWFoYXQpCiAgY29ycl9UdiA8LSBjb3JyZGV2ZGV2ZGV2ZGV2b2J2KHRfdGVzdCwgdF90ciwgZXRhaGF0MSkjIGModF90ZXN0LCB0X3RyKQogIGthcHBhX1R2IDwtIHNpZ21haGF0KmNvcnJfVHYjIGModF90ZXN0LCB0X3RyKQogIAogIGNvcnJfdG9UdiA8LSBjb3JyZGV2ZGV2ZGV2ZGV2b2J2KHRfdHIsIHRfdGVzdCwgZXRhaGF0MSkjIGModF90ciwgdF90ZXN0KQogIGthcHBhX3RvVHYgPC0gc2lnbWFoYXQgKmNvcnJfdG9UdiMgIGMoLiwgLikKICAKICBnYW1tYV92ZWMgPC0gIGthcHBhX1R2JSola2FwcGEwX2ludl9mMAogIGRrYXBwYSA8LSBzaWdtYWhhdCooMSArIGcpKnJlcCgxLCBkaW0odF90ZXN0KVsxXSpkaW0oZl90cilbMl0pCiAga2FwcGFfdmVjIDwtIChwbWF4KGRrYXBwYSAtYXBwbHkodCgoa2FwcGFfVHYpJSolKGthcHBhMF9pbnYpKSooa2FwcGFfdG9UdiksIDIsIHN1bSksIDEwKiooLTEzKSkpCiAgcmV0dXJuKGxpc3QocHJlZF9tZWFuPWdhbW1hX3ZlYywgcHJlZF92YXI9a2FwcGFfdmVjKSl9CmBgYAoKYGBge3J9CgpwcmVkbWVhbmRldmRldiA9cmVwKDAsZGltKGNlbnRlcmNvb3JkcylbMV0pCnByZWRMQmRldmRldiA9ICByZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZFVCZGV2ZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKZm9yKGsgaW4gMTpkaW0oY2VudGVyY29vcmRzKVsxXSl7CiAgcHJlZGluZm8gPC0gcHJlZGljdEdQZGV2ZGV2KGZpdHRlZF9pbmZvLHQoYXMubWF0cml4KGNlbnRlcmNvb3Jkc1trLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuZGV2ZGV2W2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQmRldmRldltrXSA9IHByZWRtZWFuZGV2ZGV2W2tdIC0gcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCiAgcHJlZFVCZGV2ZGV2W2tdID0gcHJlZG1lYW5kZXZkZXZba10rIHFub3JtKDEtMC4wNS8yKSpzcXJ0KHByZWRpbmZvJHByZWRfdmFyKQp9CnByaW50KHByZWRtZWFuZGV2ZGV2W2luZHNdKQpwcmludChwcmVkTEJkZXZkZXZbaW5kc10pCnByaW50KHByZWRVQmRldmRldltpbmRzXSkKYGBgCgoKYGBge3J9CnBsb3QoY2VudGVyY29vcmRzWywxXSxwcmVkbWVhbmRldmRldiwgeWxpbSA9IGMobWluKHByZWRMQmRldiksbWF4KHByZWRVQmRldikpKQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRMQmRldmRldikKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkVUJkZXZkZXYpCmBgYAoKYGBge3J9CgpwcmVkbWVhbmRldmRldmRldiA9cmVwKDAsZGltKGNlbnRlcmNvb3JkcylbMV0pCnByZWRMQmRldmRldmRldiA9ICByZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZFVCZGV2ZGV2ZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKZm9yKGsgaW4gMTpkaW0oY2VudGVyY29vcmRzKVsxXSl7CiAgcHJlZGluZm8gPC0gcHJlZGljdEdQZGV2ZGV2ZGV2KGZpdHRlZF9pbmZvLHQoYXMubWF0cml4KGNlbnRlcmNvb3Jkc1trLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuZGV2ZGV2ZGV2W2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQmRldmRldmRldltrXSA9IHByZWRtZWFuZGV2ZGV2ZGV2W2tdIC0gcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCiAgcHJlZFVCZGV2ZGV2ZGV2W2tdID0gcHJlZG1lYW5kZXZkZXZkZXZba10rIHFub3JtKDEtMC4wNS8yKSpzcXJ0KHByZWRpbmZvJHByZWRfdmFyKQp9CnByaW50KHByZWRtZWFuZGV2ZGV2ZGV2W2luZHNdKQpwcmludChwcmVkTEJkZXZkZXZkZXZbaW5kc10pCnByaW50KHByZWRVQmRldmRldmRldltpbmRzXSkKYGBgCgoKYGBge3J9CnBsb3QoY2VudGVyY29vcmRzWywxXSxwcmVkbWVhbmRldmRldmRldiwgeWxpbSA9IGMobWluKHByZWRMQmRldmRldmRldiksbWF4KHByZWRVQmRldmRldmRldikpKQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRMQmRldmRldmRldikKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkVUJkZXZkZXZkZXYpCmBgYAoKYGBge3J9CgpwcmVkbWVhbmRldmRldmRldmRldiA9cmVwKDAsZGltKGNlbnRlcmNvb3JkcylbMV0pCnByZWRMQmRldmRldmRldmRldiA9ICByZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZFVCZGV2ZGV2ZGV2ZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKZm9yKGsgaW4gMTpkaW0oY2VudGVyY29vcmRzKVsxXSl7CiAgcHJlZGluZm8gPC0gcHJlZGljdEdQZGV2ZGV2ZGV2ZGV2KGZpdHRlZF9pbmZvLHQoYXMubWF0cml4KGNlbnRlcmNvb3Jkc1trLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuZGV2ZGV2ZGV2ZGV2W2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQmRldmRldmRldmRldltrXSA9IHByZWRtZWFuZGV2ZGV2ZGV2ZGV2W2tdIC0gcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCiAgcHJlZFVCZGV2ZGV2ZGV2ZGV2W2tdID0gcHJlZG1lYW5kZXZkZXZkZXZkZXZba10rIHFub3JtKDEtMC4wNS8yKSpzcXJ0KHByZWRpbmZvJHByZWRfdmFyKQp9CnByaW50KHByZWRtZWFuZGV2ZGV2ZGV2ZGV2W2luZHNdKQpwcmludChwcmVkTEJkZXZkZXZkZXZkZXZbaW5kc10pCnByaW50KHByZWRVQmRldmRldmRldmRldltpbmRzXSkKYGBgCgoKYGBge3J9CnBsb3QoY2VudGVyY29vcmRzWywxXSxwcmVkbWVhbmRldmRldmRldmRldiwgeWxpbSA9IGMobWluKHByZWRMQmRldmRldmRldmRldiksbWF4KHByZWRVQmRldmRldmRldmRldikpKQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRMQmRldmRldmRldmRldikKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkVUJkZXZkZXZkZXZkZXYpCmBgYAoKYGBge3J9CmtCVCA9ICgyOTMqMC4wMDE5ODUpCnByaW50KChleHAocHJlZG1lYW5baW5kc1s0XV0tbWluKHByZWRtZWFuKSkva0JUKSpzcXJ0KC0yKnBpKmtCVC9wcmVkbWVhbmRldmRldltpbmRzWzRdXSkgKyBrQlQqKCgxLzgpKihwcmVkbWVhbmRldmRldmRldmRldltpbmRzWzRdXS8ocHJlZG1lYW5kZXZkZXZbaW5kc1s0XV0pXjIpIC0gKDUvMjQpKigocHJlZG1lYW5kZXZkZXZkZXZbaW5kc1s0XV0pXjIvKHByZWRtZWFuZGV2ZGV2W2luZHNbNF1dKV4zKSkpCnByaW50KChleHAocHJlZExCW2luZHNbNF1dLW1pbihwcmVkTEIpKS9rQlQpKnNxcnQoLTIqcGkqa0JUL3ByZWRMQmRldmRldltpbmRzWzRdXSkrIGtCVCooKDEvOCkqKHByZWRMQmRldmRldmRldmRldltpbmRzWzRdXS8ocHJlZExCZGV2ZGV2W2luZHNbNF1dKV4yKSAtICg1LzI0KSooKHByZWRMQmRldmRldmRldltpbmRzWzRdXSleMi8ocHJlZExCZGV2ZGV2W2luZHNbNF1dKV4zKSkpCgpgYGAKCmBgYHtyfQpmdHJhaW4gPC0gYXMubWF0cml4KEZyZWVFX1Byb3BhbmVkaW9sLmRmWywzXSkKaW5wdXR0cmFpbiA8LWFzLm1hdHJpeChGcmVlRV9Qcm9wYW5lZGlvbC5kZlssMToyXSkKZml0dGVkX2luZm8gPC0gZml0R1AobG9nZXRhPWMoMSwgMSwxKSwgbG93ZXJiPWMoMCwgMCwwKSwgdXBwZXJiPWMoNSwgNSw1KSwgaW5wdXR0cmFpbiwgZnRyYWluKQpmaXR0ZWRfaW5mbwpgYGAKCgpgYGB7cn0KcHJlZG1lYW4gPXJlcCgwLGRpbShpbnB1dHRyYWluKVsxXSkKcHJlZExCID0gIHJlcCgwLGRpbShpbnB1dHRyYWluKVsxXSkKcHJlZFVCID1yZXAoMCxkaW0oaW5wdXR0cmFpbilbMV0pCmZvcihrIGluIDE6ZGltKGlucHV0dHJhaW4pWzFdKXsKICBwcmVkaW5mbyA8LSBwcmVkaWN0R1AoZml0dGVkX2luZm8sdChhcy5tYXRyaXgoaW5wdXR0cmFpbltrLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuW2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQltrXSA9IHByZWRtZWFuW2tdLSBxbm9ybSgxLTAuMDUvMikqc3FydChwcmVkaW5mbyRwcmVkX3ZhcikKICBwcmVkVUJba10gPSBwcmVkbWVhbltrXSsgcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCn0KcHJpbnQocGFzdGUoJ1BDR1BNZXRob2Q6IFRlc3QgTVNFOicscm91bmQobWVhbigocHJlZG1lYW4tZnRyYWluKV4yKSwzKSwnLCBDb3ZlcmFnZTonLDEwMCpyb3VuZChtZWFuKChwcmVkTEI8ZnRyYWluKSoocHJlZFVCPmZ0cmFpbikpLDMpLCclJyxzZXAgPScnKSkKYGBgCmBgYHtyfQpGcmVlRV9Qcm9wYW5lZGlvbF9QcmVkLmRmIDwtIEZyZWVFX1Byb3BhbmVkaW9sLmRmCkZyZWVFX1Byb3BhbmVkaW9sX0xCLmRmIDwtIEZyZWVFX1Byb3BhbmVkaW9sLmRmCkZyZWVFX1Byb3BhbmVkaW9sX1VCLmRmIDwtIEZyZWVFX1Byb3BhbmVkaW9sLmRmCgoKRnJlZUVfUHJvcGFuZWRpb2xfUHJlZC5kZiRmcmVlRSA8LSBwcmVkbWVhbgpGcmVlRV9Qcm9wYW5lZGlvbF9MQi5kZiRmcmVlRSA8LSBwcmVkTEIKRnJlZUVfUHJvcGFuZWRpb2xfVUIuZGYkZnJlZUUgPC0gcHJlZFVCCgp6IDwtIGFzLm1hdHJpeChzcHJlYWQoRnJlZUVfUHJvcGFuZWRpb2wuZGYsIGBSYWRpYWwtY29vcmRgLCBmcmVlRSlbLC0xXSkKenByZWQgPC0gYXMubWF0cml4KHNwcmVhZChGcmVlRV9Qcm9wYW5lZGlvbF9QcmVkLmRmLCBgUmFkaWFsLWNvb3JkYCwgZnJlZUUpWywtMV0pCnpMQiA8LSBhcy5tYXRyaXgoc3ByZWFkKEZyZWVFX1Byb3BhbmVkaW9sX0xCLmRmLCBgUmFkaWFsLWNvb3JkYCwgZnJlZUUpWywtMV0pCnpVQiA8LSBhcy5tYXRyaXgoc3ByZWFkKEZyZWVFX1Byb3BhbmVkaW9sX1VCLmRmLCBgUmFkaWFsLWNvb3JkYCwgZnJlZUUpWywtMV0pCgpjb2xvciA8LSByZXAoMSwgbGVuZ3RoKHopKQpkaW0oY29sb3IpIDwtIGRpbSh6KQoKeCA8LXVuaXF1ZShGcmVlRV9Qcm9wYW5lZGlvbC5kZiRgUmFkaWFsLWNvb3JkYCkKeSA8LXVuaXF1ZShGcmVlRV9Qcm9wYW5lZGlvbC5kZiRgQXhpYWwtY29vcmRgKQpwIDwtIHBsb3RfbHkoKSAlPiUgYWRkX3RyYWNlKCB4PXgsIHk9eSwgej16LHR5cGUgPSAic3VyZmFjZSIsY29sb3JzY2FsZSA9IGxpc3QoYygwLCAxKSwgYygidGFuIiwgImJsdWUiKSkpCgpkZjEgPC0gc3BsaXQoRnJlZUVfUHJvcGFuZWRpb2xfTEIuZGYsIEZyZWVFX1Byb3BhbmVkaW9sX0xCLmRmJGBSYWRpYWwtY29vcmRgKQpkZjIgPC0gc3BsaXQoRnJlZUVfUHJvcGFuZWRpb2xfTEIuZGYsIEZyZWVFX1Byb3BhbmVkaW9sX0xCLmRmJGBBeGlhbC1jb29yZGApCgpmb3IoaSBpbiBzZXFfYWxvbmcoZGYxKSl7CiAgZGZfc3AgPC0gZGYxW1tpXV0KICBwIDwtIGFkZF90cmFjZShwLCBsaW5lID0gbGlzdCgKICAgIGNvbG9yID0gIiMwMDAwMDAiLCAKICAgIHdpZHRoID0gMgogICksIAogIG1vZGUgPSAibGluZXMiLCAKICB0eXBlID0gInNjYXR0ZXIzZCIsIAogIHggPSBkZl9zcCRgUmFkaWFsLWNvb3JkYCwKICB5ID0gZGZfc3AkYEF4aWFsLWNvb3JkYCwKICB6ID0gZGZfc3AkZnJlZUUsCiAgc2hvd2xlZ2VuZCA9IEZBTFNFKQp9CmZvcihpIGluIHNlcV9hbG9uZyhkZjIpKXsKICBkZl9zcCA8LSBkZjJbW2ldXQogIHAgPC0gYWRkX3RyYWNlKHAsIGxpbmUgPSBsaXN0KAogICAgY29sb3IgPSAiIzAwMDAwMCIsIAogICAgd2lkdGggPSAyCiAgKSwgCiAgbW9kZSA9ICJsaW5lcyIsIAogIHR5cGUgPSAic2NhdHRlcjNkIiwgCiAgeCA9IGRmX3NwJGBSYWRpYWwtY29vcmRgLAogIHkgPSBkZl9zcCRgQXhpYWwtY29vcmRgLAogIHogPSBkZl9zcCRmcmVlRSwKICBzaG93bGVnZW5kID0gRkFMU0UpCn0KCmRmMSA8LSBzcGxpdChGcmVlRV9Qcm9wYW5lZGlvbF9VQi5kZiwgRnJlZUVfUHJvcGFuZWRpb2xfVUIuZGYkYFJhZGlhbC1jb29yZGApCmRmMiA8LSBzcGxpdChGcmVlRV9Qcm9wYW5lZGlvbF9VQi5kZiwgRnJlZUVfUHJvcGFuZWRpb2xfVUIuZGYkYEF4aWFsLWNvb3JkYCkKCmZvcihpIGluIHNlcV9hbG9uZyhkZjEpKXsKICBkZl9zcCA8LSBkZjFbW2ldXQogIHAgPC0gYWRkX3RyYWNlKHAsIGxpbmUgPSBsaXN0KAogICAgY29sb3IgPSAiI0ZGMDAwMCIsIAogICAgd2lkdGggPSAyCiAgKSwgCiAgbW9kZSA9ICJsaW5lcyIsIAogIHR5cGUgPSAic2NhdHRlcjNkIiwgCiAgeCA9IGRmX3NwJGBSYWRpYWwtY29vcmRgLAogIHkgPSBkZl9zcCRgQXhpYWwtY29vcmRgLAogIHogPSBkZl9zcCRmcmVlRSwKICBzaG93bGVnZW5kID0gRkFMU0UpCn0KZm9yKGkgaW4gc2VxX2Fsb25nKGRmMikpewogIGRmX3NwIDwtIGRmMltbaV1dCiAgcCA8LSBhZGRfdHJhY2UocCwgbGluZSA9IGxpc3QoCiAgICBjb2xvciA9ICIjRkYwMDAwIiwgCiAgICB3aWR0aCA9IDIKICApLCAKICBtb2RlID0gImxpbmVzIiwgCiAgdHlwZSA9ICJzY2F0dGVyM2QiLCAKICB4ID0gZGZfc3AkYFJhZGlhbC1jb29yZGAsCiAgeSA9IGRmX3NwJGBBeGlhbC1jb29yZGAsCiAgeiA9IGRmX3NwJGZyZWVFLAogIHNob3dsZWdlbmQgPSBGQUxTRSkKfQpwJT4lIGxheW91dChzY2VuZSA9IGxpc3QoIHhheGlzID0gbGlzdCh0aXRsZSA9ICJSYWRpdXMsIHIiKSx5YXhpcyA9IGxpc3QodGl0bGUgPSAiQXhpYWwsIHoiKSkpCmBgYAoKYGBge3J9CgoKYXhpYWxsZW4gPC0gMTAwMApheGlhbGNvb3JkIDwtIHVuaXF1ZShGcmVlRV9Qcm9wYW5hbC5kZiRgQXhpYWwtY29vcmRgKQpjZW50ZXJjb29yZHMgPC1hcy5tYXRyaXgoIGRhdGEuZnJhbWUoc2VxKG1pbihheGlhbGNvb3JkKSxtYXgoYXhpYWxjb29yZCksbGVuZ3RoLm91dD1heGlhbGxlbikscmVwKDAsYXhpYWxsZW4pKSkKCnByZWRtZWFuID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZExCID0gIHJlcCgwLGRpbShjZW50ZXJjb29yZHMpWzFdKQpwcmVkVUIgPXJlcCgwLGRpbShjZW50ZXJjb29yZHMpWzFdKQpmb3IoayBpbiAxOmRpbShjZW50ZXJjb29yZHMpWzFdKXsKICBwcmVkaW5mbyA8LSBwcmVkaWN0R1AoZml0dGVkX2luZm8sdChhcy5tYXRyaXgoY2VudGVyY29vcmRzW2ssXSkpLCBpbnB1dHRyYWluLCBmdHJhaW4pCiAgcHJlZG1lYW5ba10gPC0gcHJlZGluZm8kcHJlZF9tZWFuCiAgcHJlZExCW2tdID0gcHJlZG1lYW5ba10tIHFub3JtKDEtMC4wNS8yKSpzcXJ0KHByZWRpbmZvJHByZWRfdmFyKQogIHByZWRVQltrXSA9IHByZWRtZWFuW2tdKyBxbm9ybSgxLTAuMDUvMikqc3FydChwcmVkaW5mbyRwcmVkX3ZhcikKfQppbmQgPSB3aGljaChwcmVkbWVhbj09bWF4KHByZWRtZWFuKSkKcHJpbnQoY2VudGVyY29vcmRzW2luZCwxXSkKcHJpbnQocHJlZG1lYW5baW5kXSkKcHJpbnQocHJlZExCW2luZF0pCnByaW50KHByZWRVQltpbmRdKQpgYGAKCmBgYHtyfQpwbG90KGNlbnRlcmNvb3Jkc1ssMV0scHJlZG1lYW4sIHlsaW0gPSBjKG1pbihwcmVkTEIpLG1heChwcmVkVUIpKSkKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkTEIpCmxpbmVzKGNlbnRlcmNvb3Jkc1ssMV0scHJlZFVCKQoKYGBgCgpgYGB7cn0KCnByZWRtZWFuZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZExCZGV2ID0gIHJlcCgwLGRpbShjZW50ZXJjb29yZHMpWzFdKQpwcmVkVUJkZXYgPXJlcCgwLGRpbShjZW50ZXJjb29yZHMpWzFdKQpmb3IoayBpbiAxOmRpbShjZW50ZXJjb29yZHMpWzFdKXsKICBwcmVkaW5mbyA8LSBwcmVkaWN0R1BkZXYoZml0dGVkX2luZm8sdChhcy5tYXRyaXgoY2VudGVyY29vcmRzW2ssXSkpLCBpbnB1dHRyYWluLCBmdHJhaW4pCiAgcHJlZG1lYW5kZXZba10gPC0gcHJlZGluZm8kcHJlZF9tZWFuCiAgcHJlZExCZGV2W2tdID0gcHJlZG1lYW5kZXZba10gLSBxbm9ybSgxLTAuMDUvMikqc3FydChwcmVkaW5mbyRwcmVkX3ZhcikKICBwcmVkVUJkZXZba10gPSBwcmVkbWVhbmRldltrXSsgcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCn0KaW5kcyA9IHdoaWNoKGRpZmYoc2lnbihwcmVkbWVhbmRldikpICE9IDApCnByaW50KGNlbnRlcmNvb3Jkc1tpbmRzLDFdKQoKcHJpbnQocHJlZG1lYW5kZXZbaW5kc10pCnByaW50KHByZWRMQmRldltpbmRzXSkKcHJpbnQocHJlZFVCZGV2W2luZHNdKQpgYGAKCgpgYGB7cn0KcGxvdChjZW50ZXJjb29yZHNbLDFdLHByZWRtZWFuZGV2LCB5bGltID0gYyhtaW4ocHJlZExCZGV2KSxtYXgocHJlZFVCZGV2KSkpCmxpbmVzKGNlbnRlcmNvb3Jkc1ssMV0scHJlZExCZGV2KQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRVQmRldikKCmBgYAoKYGBge3J9CgpwcmVkbWVhbmRldmRldiA9cmVwKDAsZGltKGNlbnRlcmNvb3JkcylbMV0pCnByZWRMQmRldmRldiA9ICByZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZFVCZGV2ZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKZm9yKGsgaW4gMTpkaW0oY2VudGVyY29vcmRzKVsxXSl7CiAgcHJlZGluZm8gPC0gcHJlZGljdEdQZGV2ZGV2KGZpdHRlZF9pbmZvLHQoYXMubWF0cml4KGNlbnRlcmNvb3Jkc1trLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuZGV2ZGV2W2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQmRldmRldltrXSA9IHByZWRtZWFuZGV2ZGV2W2tdIC0gcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCiAgcHJlZFVCZGV2ZGV2W2tdID0gcHJlZG1lYW5kZXZkZXZba10rIHFub3JtKDEtMC4wNS8yKSpzcXJ0KHByZWRpbmZvJHByZWRfdmFyKQp9CnByaW50KHByZWRtZWFuZGV2ZGV2W2luZHNdKQpwcmludChwcmVkTEJkZXZkZXZbaW5kc10pCnByaW50KHByZWRVQmRldmRldltpbmRzXSkKYGBgCgoKYGBge3J9CnBsb3QoY2VudGVyY29vcmRzWywxXSxwcmVkbWVhbmRldmRldiwgeWxpbSA9IGMobWluKHByZWRMQmRldiksbWF4KHByZWRVQmRldikpKQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRMQmRldmRldikKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkVUJkZXZkZXYpCmBgYAoKYGBge3J9CgpwcmVkbWVhbmRldmRldmRldiA9cmVwKDAsZGltKGNlbnRlcmNvb3JkcylbMV0pCnByZWRMQmRldmRldmRldiA9ICByZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZFVCZGV2ZGV2ZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKZm9yKGsgaW4gMTpkaW0oY2VudGVyY29vcmRzKVsxXSl7CiAgcHJlZGluZm8gPC0gcHJlZGljdEdQZGV2ZGV2ZGV2KGZpdHRlZF9pbmZvLHQoYXMubWF0cml4KGNlbnRlcmNvb3Jkc1trLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuZGV2ZGV2ZGV2W2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQmRldmRldmRldltrXSA9IHByZWRtZWFuZGV2ZGV2ZGV2W2tdIC0gcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCiAgcHJlZFVCZGV2ZGV2ZGV2W2tdID0gcHJlZG1lYW5kZXZkZXZkZXZba10rIHFub3JtKDEtMC4wNS8yKSpzcXJ0KHByZWRpbmZvJHByZWRfdmFyKQp9CnByaW50KHByZWRtZWFuZGV2ZGV2ZGV2W2luZHNdKQpwcmludChwcmVkTEJkZXZkZXZkZXZbaW5kc10pCnByaW50KHByZWRVQmRldmRldmRldltpbmRzXSkKYGBgCgoKYGBge3J9CnBsb3QoY2VudGVyY29vcmRzWywxXSxwcmVkbWVhbmRldmRldmRldiwgeWxpbSA9IGMobWluKHByZWRMQmRldmRldmRldiksbWF4KHByZWRVQmRldmRldmRldikpKQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRMQmRldmRldmRldikKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkVUJkZXZkZXZkZXYpCmBgYAoKYGBge3J9CgpwcmVkbWVhbmRldmRldmRldmRldiA9cmVwKDAsZGltKGNlbnRlcmNvb3JkcylbMV0pCnByZWRMQmRldmRldmRldmRldiA9ICByZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKcHJlZFVCZGV2ZGV2ZGV2ZGV2ID1yZXAoMCxkaW0oY2VudGVyY29vcmRzKVsxXSkKZm9yKGsgaW4gMTpkaW0oY2VudGVyY29vcmRzKVsxXSl7CiAgcHJlZGluZm8gPC0gcHJlZGljdEdQZGV2ZGV2ZGV2ZGV2KGZpdHRlZF9pbmZvLHQoYXMubWF0cml4KGNlbnRlcmNvb3Jkc1trLF0pKSwgaW5wdXR0cmFpbiwgZnRyYWluKQogIHByZWRtZWFuZGV2ZGV2ZGV2ZGV2W2tdIDwtIHByZWRpbmZvJHByZWRfbWVhbgogIHByZWRMQmRldmRldmRldmRldltrXSA9IHByZWRtZWFuZGV2ZGV2ZGV2ZGV2W2tdIC0gcW5vcm0oMS0wLjA1LzIpKnNxcnQocHJlZGluZm8kcHJlZF92YXIpCiAgcHJlZFVCZGV2ZGV2ZGV2ZGV2W2tdID0gcHJlZG1lYW5kZXZkZXZkZXZkZXZba10rIHFub3JtKDEtMC4wNS8yKSpzcXJ0KHByZWRpbmZvJHByZWRfdmFyKQp9CnByaW50KHByZWRtZWFuZGV2ZGV2ZGV2ZGV2W2luZHNdKQpwcmludChwcmVkTEJkZXZkZXZkZXZkZXZbaW5kc10pCnByaW50KHByZWRVQmRldmRldmRldmRldltpbmRzXSkKYGBgCgoKYGBge3J9CnBsb3QoY2VudGVyY29vcmRzWywxXSxwcmVkbWVhbmRldmRldmRldmRldiwgeWxpbSA9IGMobWluKHByZWRMQmRldmRldmRldmRldiksbWF4KHByZWRVQmRldmRldmRldmRldikpKQpsaW5lcyhjZW50ZXJjb29yZHNbLDFdLHByZWRMQmRldmRldmRldmRldikKbGluZXMoY2VudGVyY29vcmRzWywxXSxwcmVkVUJkZXZkZXZkZXZkZXYpCmBgYAoKYGBge3J9CmtCVCA9ICgyOTMqMC4wMDE5ODUpCnByaW50KChleHAocHJlZG1lYW5baW5kc1s0XV0tbWluKHByZWRtZWFuKSkva0JUKSpzcXJ0KC0yKnBpKmtCVC9wcmVkbWVhbmRldmRldltpbmRzWzRdXSkgKyBrQlQqKCgxLzgpKihwcmVkbWVhbmRldmRldmRldmRldltpbmRzWzRdXS8ocHJlZG1lYW5kZXZkZXZbaW5kc1s0XV0pXjIpIC0gKDUvMjQpKigocHJlZG1lYW5kZXZkZXZkZXZbaW5kc1s0XV0pXjIvKHByZWRtZWFuZGV2ZGV2W2luZHNbNF1dKV4zKSkpCnByaW50KChleHAocHJlZExCW2luZHNbNF1dLW1pbihwcmVkTEIpKS9rQlQpKnNxcnQoLTIqcGkqa0JUL3ByZWRMQmRldmRldltpbmRzWzRdXSkrIGtCVCooKDEvOCkqKHByZWRMQmRldmRldmRldmRldltpbmRzWzRdXS8ocHJlZExCZGV2ZGV2W2luZHNbNF1dKV4yKSAtICg1LzI0KSooKHByZWRMQmRldmRldmRldltpbmRzWzRdXSleMi8ocHJlZExCZGV2ZGV2W2luZHNbNF1dKV4zKSkpCgpgYGA=